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ABSTRACT 



An application of two-dimensional Fourier Transform correlation 
for the attitude measurement of three-dimensional objects is considered. 
Included are geometric optical equations and geometric translation and 
rotation equations necessary for the image error calculations. Mini- 
mum attitude resolution equations are developed to estimate image 
quantization necessary for a given measurement error. These techniques 
are then applied as an additional measurement to the radar tracking 
problem. 
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I. INTRODUCTION 



The development of sophisticated naval weapons has resulted in the 
requirement for high accuracy measurement tests and discrete analysis of 
complex system parameters. Typical of this precise test and evaluation 
requirement is attitude measurement of three-dimensional airborne objects 
by processing multiple two-dimensional images. 

The most common method used to estimate attitude is by making 
absolute or precisely relative three-dimensional geometric trajectory 
measurements of target object extreme structural features (nose, tail, 
right wing tip, left wing tip, etc.) from camera images. By appropriate 
camera, range coordinate system, and target coordinate system trans- 
formations these measurements will yield an attitude measurement. 

A more direct method of attitude measurement involves the correlation, 
or comparison, of an image of an object at an unknown attitude with a 
second image of a similar object at a known attitude. At the Naval 
Weapons Center, an operator controlled film and television correlation 
assessment technique (FATCAT) facility use this more direct technique 
for attitude estimation [1]. 

This thesis presents an enhancement of the correlation/comparison 
method of attitude measurement by fully automating the technique incorp- 
orating a pattern recognition method of the directly digitized image. 

Image theory, resolution, and filtering leading to the two dimensional 
matched filter attitude measurement concept is discussed and this concept 
is applied as an additional measurement to a radar tracking problem. 
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II. IMAGE THEORY 



Three-dimensional to two-dimensional projections are obtained by 
viewing the real world as recorded by a camera or other imaging device. 
Image theory defines projection as the transformation of each point in 
the field of view, through a focal point, and onto a plane. This is 
also equivalent to locating the focal plane in front of the focal point 
and changing the sign of the projected points [2]. These geometrical 
relationships are illustrated in Fig. 1. The transformation of point 
X, y, z in the X, Y, Z coordinate system of Fig. 1 is given by the 
following equations: 



u = 



y + f 



( 1 ) 



V = 



y + f 



( 2 ) 



Where u = projected x dimension 
V = projected z dimension 

If f << y, then the following approximations hold; 



fx 



u = 



= kx 



(3) 



fz 

V = — ky 



(4) 
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Therefore, for a given focal length f, and target range y, a constant 
image scaling factor k, is calculated. 




A second coordinate system X', Y', V can be used to describe a 
three-dimensional object in the viewing space. If the origins of these 
two coordinate systems are coincident, and if all the axes are parallel, 
then the transformation of points from one coordinate system to the 
other would be simple. 
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If the axes of the two coordinate systems are rotated by an elevation 
angle Q around the camera Y axis, by an azimuth angle around the camera 
Z axis, and by roll angle <t> around the camera X axis, the transformation 
equation would be [3], [4]: 



’x'" 
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( 6 ) 



If in addition there was a translation of (A,B,C) of the origin of the 
object coordinate system relative to the origin of the camera coordinate 
system, the transformation equation would be: 




where 
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III. IMAGE RESOLUTION 



The three-dimensional object can be considered to be an aircraft, 
with maximum orthogonal dimensions on the focal plane of Xmax, Ymax, 

Zmax. These dimensions are illustrated on the top and side views of 
Fig. 2. Since the image is to be digitaized, the discrete maximum 
dimensions can be given in terms of pixel sampling intervals, x and y. 

A pixel is defined as the smallest unit sampling area (x,y) of resolution 
for the digitized image. 



Xmax (quantized) = XN 


( 8 ) 


Ymax(quantized) = YN 


(9) 


Zmax (quantized) = ZN 


( 10 ) 



These dimensions are illustrated in Fig, 3. The quantization of the 
image would have a minimum rotation that could be detected by a 1 pixel 
change at extreme points in the image and would be equivalent to the 
following angle: 



tan 0 



J_ 

KN 



01 ) 



This relationship is illustrated in Fig. 4. KN would be equal to XN, YN, 
or ZN depending on the angle of rotation. 
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Figure 2 

Maximum Image Dimensions, Orthogonal Views 
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□ = 1 pixel 





Figure 3 

Discrete, Orthogonal Views 
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If the axis of rotation was not orthogonal to the focal plane, the 
minimum detectable rotation of an equal sized object would be much 
larger, as indicated in the top view of Fig. 5. 



cos 0 



min 



KN-1 

KN 



KN 



1 



1-cos 0 



min 



( 12 ) 

(13) 



Sample values of KN from equations (11) and (13) for minimum detectable 
rotation are given in Table 1. 
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Side View Geometry 




Figure 5. Minimum Detectable Rotation 

Rotation Axis in Viewing Plane 
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TABLE 1 



Axis Orthogonal Axis Parallel 

to image plane (11) to image plane (13) 
0 (degrees) KN (pixels) KN (pixels) 



1 


58 


6566 


5 


12 


263 


10 


6 


66 



If it is assumed that the rotational angles are measured by lines 
through the two extreme pixels, the geometry and error sensitivity would 
be as shown in Fig. 6. For a maximum length of KN. and a 1 pixel rota- 
tion: 



1 

tan 9 = (11 ) 

KN 



For horizontal errors H: 



1 

0 = arctan z (14) 

KN - H 



1 

de = d arctan z (15) 

KN - H 



For a change of variables: 



1 

m = (16) 

KN +H 
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dH 



(17) 



dm 



(KN t H)2 



1 

de = dm 

1 + m^ 

dH 



(KN - H)^+l 



(18) 

(19) 



For example, if KN = 58, H = 1, then de/dH = 0.0003077 radians/pixel = 
0.01763 degrees/pixel. This sensitivity for horizontal errors is very 
small, as expected. 

For vertical errors: 



9 = arctan 




KN 




KN 



dV 

dm = 

KN 



dV KN 

KN^ + (1 - V)^ 



( 20 ) 



( 21 ) 

( 22 ) 



(23) 



For example, if V = 1 pixel, and KN = 58 pixels, de/dV = 0.01722 radians/ 
pixel = 0.99 degrees/pixel. This is much larger than the horizontal 
error, as might be expected by observing the geometry of Fig. 6. 
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KN pixels 



— = Reference Axis 
H = Horizontal Error 
V = Vertical Error 






<■ 




Figure 6. 

One-Dimensional Rotation Errors 
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This section has derived equations (11) and (12) to predict the 
minimum detectable rotation of a three-dimensional object. These 
equations apply to a three-dimensional object with a two-dimensional 
maximum projected dimension of KN pixels. Equation (11) applies if 
the axis of rotation of the three-dimensional object is parallel 
with the camera axis. Equation (13) applies if the axis of rotation 
of the three-dimensional object is orthogonal to the camera axis. 
Equations (19) and (23) predict the rotational errors resulting from 
horizontal and vertical measurement errors in extreme points of the 
two-dimensional image. 
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IV. IMAGE FILTERING 



Two-dimensional matched filtering can be considered to be an 
extension of one-dimensional correlation filtering [5]. The out- 
put of the two-dimensional matched filter will not be a transformed 
image with some type of improvement or enhancement, but a correlation 
plane whose amplitude at a given location corresponds to the degree 
of correlation of the input image (signal + noise) with the desired 
image (signal). This is equivalent to sliding one photographic 
transparency across another and adding up each black point that was 
aligned with another black point on the second image, each white 
point that was aligned with another white point on the second image, 
and each grey point that was aligned with another grey point of 
equal shade on the second image. 

Several practical factors must be considered in the actual 
simulation of the matched filter. The input image and the refer- 
ence or stored image could both be compared (correlated) as contin- 
uous grey scale images, but this would be a tremendous storage and 
computational task. Therefore, the reference image in the computer 
simulation is digitized and stored as discrete X,Y,Z points (FACIT 
Model) in a three-dimensional coordinate system [6]. It is recon- 
structed, after rotation and translation, as unit amplitude lines 
projected on the focal plane. It is also assumed that f << y, and 
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f and y of equations (1) and (2) are known. This knowledge of 
focal length and target range yields a constant image scaling fac- 
tor (k). This type of projection is illustrated in Fig. 7. An 
example of an aircraft FACIT Model is given in Fig. 8. 

The input grey-level image is also transformed into a similar 
format by use of an edge detector to extract outline information 
[7].' The edge of the image is defined as the boundary between the 
two regions of different grey level. 

The reference image and the transformed input image are applied 
to a matched filter to determine the degree of correlation. An 
ideal correlation function C(x,y,k,0,t,<l>) would have an impulse 
response that could be modeled by delta, (sin x/x), triangle or 
exponential functions. The impulse would occur when the reference 
and transformed images were: 

1. aligned on the image plane (x,y), 

2. equally scaled (k) and 

3. at the same attitude ,<()). 

However, due to the presence of noise and the approximation of the 
three-dimensional object with a finite element model, the correla- 
tion response is not delta, (sin x)/x, triangle or exponential in 
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Figure 8. FACIT Model of F-102 
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shape. It is not necessary to know the exact correlation response 
function or even the constant component due to a nonzero mean value 
of the reference and input image variable function. It is only 
necessary to know the values of the variables x, y, k, e , and <(> 
at peak correlation. 

The following simplifying assumptions are made for an attitude 
only correlation: 

1. the images have been scaled to an equal size (k) 

2. the peak correlation on the image plane will indicate 
translational alignment (x,y) 

Fig. 9 illustrates a typical x-y correlation between a noisy image 
and computer-generated reference image. The peak correlation at 
(0.0, 0.0) pixels indicates that both images are aligned in the 
x-y plane. This is because the noisy image was generated from a 
computer reference image. 

The attitude of the input image is unknown. Therefore, the 
measurement of the attitude of the input image is equivalent to 
applying the transformed image to a bank of matched filters with each 
"tuned" to a specific attitude. This process is illustrated in 
Fig. 10. 
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X-Y CORRELATION PLANE 




Sample 64 X 64 Pixel X-Y Correlation Plane Obtained 
from 32 X 32 Pixel Image and 32 X 32 Pixel Reference 
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max 



•^max 



•^max 



Figure 10. Image Processing Flov/ Diagram 



If there was a priori knowledge of azimuth 6 , elevation ^ , and 
roll (p , to within 10 degrees each, and a 1 degree resolution of 
measurement was required, a brute force approach would require 
10 X 10 X 10 = 1000 correlations or 1000 matched filters. Three axis 
quadratic interpolation is used to reduce the number of correlations 
required for the given resolution of measurement [3]. 

It is assumed that the variables azimuth, elevation, and roll 
are uncoupled in the joint correlation function C(x,y,6,ik <P). Assum- 
ing that there is a near orthogonal relationship allows C(x,y,e ,({> ) 
to be sequentially maximized (minimized) with the independent vari- 
ables: X, y, azimuth, elevation, and roll. 
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The pitch correlation function is modeled as a quadratic func- 
tion C(9) : 



C(e) =ae +be+c 



dC 

de 



= 2 a 0 + b 



(24) 

(25) 



Equation (25) is set to zero to maximize (minimize) C(o). 



max 



2 a 



(26) 



01 = 0 degrees, 02 = 0i + A0 degrees, and 03 = 0i - A0 degrees are 
substituted in equation (24) to solve for the unknowns a, b, and c. 
The solution of equation (26) is: 



Q = (C(0,) - C(0,)) A0 

max 2 C( 02 ) + 20 ( 03 ) - 4C(0i) 



(27) 



Similar equations are written for C(ii»), and C(ij)). 

The computer simulations listed in Appendix B implement equa- 
ion (27). The three correlations C( 0 i), 6 ( 02 ) and 6 ( 03 ) are calculated 
with A0 = 5 degrees for a rough approximation of Then three 
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correlations C( 6 i ), C (02 ) and C (93 ) are calculated with A9 = 2 
degrees for a sharper quadratic approximation of . In addition, 
the and degree intervals are tested to make sure that the mid- 
point correlation C( 0 i) was greater than or equal to the end point 
correlations 0 ( 62 ) and 0 ( 63 ). If either is not, the interval is shifted 
by 2 or 5 degrees in the direction of the larger correlation and the 
mid-point is tested again. This dual pass computation reduces the 
unwanted coupling effect of the attitude variables. The +5 degree 
calculation detects the general angle of maximum correlation. The 
j ^2 degree calculation reduces the estimated angle error by not being 
as sensitive to unsymmetry and nonlinearities in the correlation func- 
tion. 



A sample of the dual pass azimuth estimation is illustrated in 
Fig. 11. For this sample, noisy edge data of a simulated aircraft at 
an azimuth angle of -45.0 degrees is used as input data. This data 
is correlated with computer-generated edge data from a stored air- 
craft model. The initial estimate of attitude is -48.0 degrees. The 
final azimuth calculated is -44.11 degrees, a 0.89 degree error from 
the true -45.0 degree aircraft azimuth attitude. This is a typical 
maximum error. For example, if the initial estimate of azimuth is 
-42.0 degrees, the final error would be 0.0 degrees. 
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LEGEND 



□ = CORRELATION DATA 
o = FIRST INTERPOLATION 




Figure 11 

Azimuth Correlation Function 
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This section has introduced the use of two-dimensional correla- 
tion to determine the attitude of a three-dimensional object. Equation 
(24) is used as a second order polynomial model for the correlation 
function of each attitude variable. This approximation reduces the 
number of computations necessary to determine the peak correlation for 
each attitude variable. This quadratic approximation was demonstrated 
with a computer simulation that tested for the maximum correlation of 
an image as a function of attitude. The simulated image was previously 
generated by adding edge noise to a computer-generated image at a known 
attitude of -45.0 degrees. By starting the correlation algorithm with 
a -3.0 degree offset, the final error was found to be 0.89 degrees. 
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V, MATCHED FILTER EQUATION 



It is well known that one-dimension correlation can be calculated 
by a simple multiplication of two functions that are Fourier trans- 
formed [ 9 ]. 

” .Z ^ h(t+T) dt = x(t) * h(-t) (28) 

where t is the dummy variable of integration and * is used to 
symbolize the convolution integral. 

C(w) = X(w) H*(w) = X*(w) H(w) (29) 

where w = 2Trf = frequency in radians and * is used to indicate com- 
plex conjugation; that is, if 

H(w) = R(w) + jl(w), then H*(w) = R(w) - jl(w). 

c(t) <=> C(w) (30) 

where <=> is used to indicate a Fourier transform pair; that is 

C(w) is the Fourier transform of c(t). 
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This same relationship can be shown in two-dimensions and in discrete 
form [5]. 



N-1 M-1 

c(k»l) = s z X (i,j) h(k+i,l+j) (31) 

i=0 j=0 



C(m,n) = X*(m,n) H(m,n) 



(32) 



C(m,n) = 



K-1 L-1 

,• 2 nmk i 2 nnl 

z xCk.De '-* K e L 

k=0 1=0 



I-l J-1 

. 2 II mi 

z- z h(i,j)e ^ I 
1=0 j =0 



e 



. 2 n nj 

3 ~ 1 } 



(33) 



m = variable in first dimension 



n = variable in second dimension 



K = I = M = Period in first dimension 



1 _ = j = N = Period in second dimension 
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A complete listing of the transformation indicated in equations (32) 
and (33) is given in Appendix A. 

The reason the Fourier domain correlation approach is used is 
that for large data arrays the Fast Fourier Transform (FFT) computa- 
tion time is less than that of direct time or space domain correlation. 
Also, Fourier transform and array processing hardware is available for 
many scientific computers. 

The continuous correlation equation (28) states that an infinite 
integration is necessary to determine the correlation at a particular 
time. This is not the case for finite sampled functions. The summa- 
tion interval for the sampled functions requires that the majority of 
the information in the two functions be contained in the finite sampled 
interval for the correlation to approximate the continuous correlation. 
When x(t) and h(t) are sampled with an impulse function A(t) over a 
finite interval, 0<t<t„^^ = NT, their Fourier transforms x(w)*a(w) 
and H(w)*a(w) will be periodic with period w = 1/T, where T = sampling 
interval of the impulse function A(t). When the continuous function 
x(t) is sampled over the interval 0< and Fourier transformed, 

this is equivalent to calculating the Fourier Series coefficients of a 
sampled periodic sequence x (nT) of period N, where N = t Jl and 

r r 

T = time interval between samples. If Xp(nT) is time shifted in 
relation to x(nT) by m samples, and the first N samples are trans- 
formed, it can be shown that the resulting complex Fourier transform 
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win change in phase but not in amplitude or period [9]- In the 
time domain, this shift of m samples is due to the periodicity 
of Xp (nT) and is referred to as a circular shift of a sequence [5]. 

The sequence appears to be circular because the last sample of x(nT) 
that is shifted beyond the interval N T appears at the beginning of 
the sequence. Likewise, in the two-dimensional correlation there is a 
circular shift in two dimensions when one image x(i,j) is correlated 
with another h(i,j). The two-dimensional periodic or circular sampled 
function is equivalent to joining parallel horizontal and vertical edges 
of each image and rotating one image on each of its axes while the 
other image remains stationary. This two-dimensional circular shifting 
is depicted in Fig. 12, In this figure, image 2 has shifted n samples 
in the horizontal direction with respect to the fixed image 1. Image 2 
has shifted m samples in the vertical direction with respect to the 
fixed image 1, To clarify that image 2 has circularly "wrapped around" 
image 1, the corners are labeled A, B, C and D, 
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0 



The continuous images 1 and 2 are sampled at N and M discrete 
points in the x and y axis respectively. This sampling results in 
a periodic equation (33) in the transformed frequency domain of per- 
iod N and M. If the correlation product of the two sampled images 
has nonzero terms beyond the N and M dimensions, there will be an 
overlap of nonzero samples in the frequency domain. This overlap 
phenomenon of a high frequency component taking on the identity of 
a lower frequency component is referred to as frequency aliasing. 
Likewise, spatial aliasing would occur if there was an overlap of the 
periodic spatial functions; for example, if the first and last samples 
added. 

To test the correlation method of attitude measurement, a com- 
puter simulation program, FAC listed in Appendix B, was run. Discrete 
boundary points were generated as input data for a known aircraft 
attitude as illustrated in Fig. 13. A noisy test object was generated 
that was within +^5 degrees in azimuth, elevation and roll of the ref- 
erence object. The noisy boundary, as described in the next paragraph, 
was compared with the noise-free computer-generated edge data to esti- 
mate the aircraft attitude. The representative errors of attitude 
(azimuth, elevation and roll) were computed by comparing the attitude 
that was iteratively predicted by equation (27) with the actual atti- 
tude. These results are summarized in Table II. 
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Figure 13 

Boundary Points of FACIT Model 
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Table II. 



Representative Attitude Errors (Degrees) 

Three Sample Average, 1 Pixel Noise Standard Deviation 



Azimuth 


0.12 


Elevation 


0.55 


Roll 


0.71 



It was noted that the errors increased with each attitude measurement. 
This is because of the sequential coupling of a previous attitude 
angle error to the next quadratic angle interpolation. 

The quadratic correlation attitude estimator (27) was initial- 
ized with data to within 10 degrees of the known aircraft attitude. 
Boundary points on a 32 x 32 grid were generated from the two- 
dimensional projection of the test aircraft. The number of boundary 
points generated were a function of the test aircraft attitude, but 
typically 108 points were generated: 

(27 each: maximum horizontal, minimum horizontal, maximum 

vertical, minimum vertical). 
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To simulate detector noise, Gaussian noise with zero mean and one 
pixel standard deviation was added to each boundary point, horizon- 
tally and vertically. Since the boundary point function was on a 
discrete 32 x 32 grid, this had the effect of quantizing the noise 
to unit pixel intervals. For example, if the absolute value of the 
Gaussian noise was less than 0.5 pixel, no noise would be added. If 
the absolute value of the Gaussian noise was between 0.5 and 1.5 
pixels, one pixel of noise was added. Likewise, two or more pixels 
of noise would be added for greater outputs of the generator. These 
noisy boundary points are shown in Fig. 14. 

An uncalibrated system test with real data was made from a 
digitized video image of a F102 aircraft. The intensity data from 
the image was transformed into edge points. The edge points were 
correlated with the model data stored in the computer memory. This 
predicted the aircraft's attitude. 

The x,y resolution of the digitization was 340 pixel columns 
X 220 pixel rows and the intensity resolution was 1 to 256. As the 
photograph in Fig. 15 shows, the aircraft image makes up only a small 
portion of the video image. Therefore, a 64 x 64 pixel window shown 
in Fig. 16 was contrast enhanced, in that the original intensity 
range of 103 to 156 was linearly transformed to a 1 to 256 intensity 
level. It is noted that the digitized image is not displayed at the 
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Figure 15 

Photograph of FI 02 Aircraft 
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Figure 16 

Digitized 64 X 64 Pixel Image of FI 02 Aircraft 
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correct aspect ratio due to the unequal sampling ratios taken verti- 
cally and horizontally. 

The edge data from Rosenfeld's Detector [7] was applied to a 
variable threshold iteration that was adjusted until greater than 
200 pixels were above the threshold. This adjustment was necessary 
so that low contrast edge pixels typically at wing, nose and tail 
extremities would exceed the threshold. These extremities are most 
important in obtaining a given attitude accuracy. Fig. 17 shows the 
results of the best computer-generated match of the stored image 
superimposed on the edge data obtained from the video image. The 
computer predicted the attitude to be; 

Azimuth = -115.0 degrees 

Elevation = 45.0 degrees 

Roll = - 5.0 degrees 

There was no other data available to establish the attitude of the 
FI 02 aircraft at the instant the camera recorded the image. It is 
presumed that further investigation would show this attitude to be 
correct to within the limits of the image resolution and the edge 
noise present. 
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Figure 17 

Best Fit of Stored Image Pixels to Video Outline Pixels of Fig. 16 
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The lack of sharpness in the transformed edge data shown in 
Fig. 17 can be attributed to several factors; 

1. Low x,y resolution (64 x 64), 

2. Low intensity range (103 - 156), 

3. Uneven lighting and shading on edges. 

4. A simple edge detector. 

The last factor would be the only one that could be significantly 
changed given a free-flying aircraft, fixed camera locations, and 
natural lighting. A review of methods used for extracting features 
based on edges and contours is included in a survey by Lavine [11]. 
More recent theories and experiments on edge detectors have been 
demonstrated by Abdou and Pratt [12]. 

This section has introduced the use of two-dimensional Fourier 
domain correlation as described by equation (32). Representative 
results of three computer correlation simulations were listed. The 
simulation compared a computer generated edge model to a previous 
computer generated model with added edge noise. Finally, a sample 
correlation of the computer generated edge model with the edge data 
detected from a real video image of an aircraft at an unknown atti- 
tude was described. 
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VI. AIDED AIRCRAFT TRACKING 



An example of the utility of the attitude measurement is demon- 
strated by utilizing an aircraft positional estimating Kalman Filter 
with combined radar and electro-optic sensor measurements [10]. It 
is difficult for the radar-only estimator to keep a close target in 
track if it is highly maneuverable. The criterion of "highly man- 
euverable" for a particular radar is dependent on the following: 

1. The sampling interval of the radar. 

2. The range gate period, or the sampling aperture. 

3. The receiver sensitivity. 

4. The return signal to noise ratio. 

5. The mount azimuth and elevation dynamics, or time constants. 

6. The measurement accuracy of azimuth, elevation and range. 

7. The point source tracking of a finite length target. 

These parameters and others contribute to the ability of the radar 
to track with minimum error. By modeling the aircraft target to be 
driven by an acceleration function that is dependent upon aircraft 
attitude, more accurate range, range rate and range acceleration state 
estimates are possible. The acceleration function is based on the 
well-known aerodynamic coupling between the target attitude and the 
direction of the target acceleration. This coupling is summarized 
in the following comments: 
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1 . 



Major accelerations (lift) are normal to the velocity vector 
and the pitch axis. 

2. Positive lift is more likely than negative lift due to pilot 

physiological factors and to structural loading design. 

3. Accelerations in the velocity direction (drag/thrust) are 

generally smaller in magnitude and of shorter duration than 
the lift accelerations. 

A complete multiple measurement system using the previously discussed 
imaging attitude measurement concept as a radar tracking aid to the 
positional estimating Kalman Filter is shown in Fig. 18. In the upper 
left corner of Fig. 18 there are n radars that provide range, azimuth 
and elevation measurements to the Iterated Extended Kalman target 
estimator. In the lower left corner there are m video cameras that 
input data to m correlators. These correlators first scale the stored 
model to a size equal to the predicted input image size. The scaling 
is a function of target. range and system gains such as telescope mag- 
nification and image pixel sampling area. The correlators then apply 
equation (33) as described in the previous section for the computer 
simulation FAC. The attitude measurements azimuth, elevation, and 
roll, from the correlators, are transformed into a common coordinate 
system. The predicted attitude of a model aircraft in a coordinated 
turn is combined with the m attitude measurements in a Kalman Filter. 
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Figure 18. Attitude Measurement and Tracking System 
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The estimated attitude from the Kalman Filter predicts an accelera- 
tion vector that is used as an input to the target estimator. 

The simulation of a target range (R) estimation was performed 
with the program KAL listed in Appendix B. The simulation KAL was a 
limited example of the radar tracking filter as illustrated in Fig. 18. 
This simplification was made possible by not including some tracking 
variables, such as azimuth and elevation, to the full three-dimensional 
trajectory filter. The trajectory was truly three-dimensional, hence 
there were changes in azimuth and elevation, but they were not used 
as system measurements or estimated as parameters. 

The Iterated Extended Kalman Filtering subroutine FTKALI is a 
modification of the IMSL Kalman Filtering subroutine FKALM [14]. 

The system was modeled as a radar track of a nominal constant accel- 
eration target. The assumptions for this trajectory are; 

1. The average acceleration rate is zero. 

2. The acceleration rate between radar scans is constant. 

3. The acceleration between scanning intervals in uncorrelated. 

The above assumptions are described in the well-known Alpha-Beta- 
Gamma Filter equations [13]- Written in a discrete form, these 
equations are: 
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( 34 ) 




(35) 




(36) 



where ar is a random sequence with zero-mean and a mean-square 
acceleration rate of 



The target trajectory was simulated in three parts: 

1. A target constant closing velocity of 100 or 175 
yards/second at the same altitude as the radar. 

2. A 10-g target turn in the azimuth plane of the radar. 

3. A 10-g target turn "up" in the geometric plane perpendicular 




ar 



2 

0 for j = k 



(37) 



0 otherwise 



to the azimuth plane and the ground projection of the radar 



range vector. 
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The final, or third turn, starts when the target has rotated through 
90.0 degrees of the 10-g turn in the azimuth plane. The maneuver 
was started at three initial ranges: 



Inputs to the filter of two noisy radar range measurements and an 
optional radial component of the predicted normal acceleration vector 
were generated at a sample rate of 10/second. A one sigma standard 
deviation Gaussian noise of 1.0 yard was added to the two range 
measurements. These inputs were input to the range filter. 

Typical results of range error for one simulation of all man- 
euvers are listed in Table III. The addition of the predicted normal 
.acceleration measurement to the radar range measurements reduced the 
peak and RMS errors for all maneuvers. The reduction in error was 
most dramatic for the close-in maneuver because of the range reversal 
for the particular azimuth and elevation geometries. At the mid- 
range and distant range trajectories there was no range reversal 
so the additional acceleration measurement was less helpful in cor- 
recting the range and range rate estimates In addition, the high 
radar sample rates and the averaging effect of using two radar measure- 
ments reduced the difference in error between the radar-only estimate 
and the combined estimate. 



1. Distant 



11275.0 yards 
1000.0 yards 
150.0 yards 



2. Mid-range 



3. Close-in 
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Table III. 



Range Error (RMS Average/Peak - in Yards) 



Maneuver 


Combined Measurement 
(Normal Acceleration 
and Two Radars) 


Two Radars Only 


1. Distant 


0.76/3.22 


0.99/4.82 


2. Mid-Range 


0.89/4.25 


1.06/5.53 


3. Close-In 


4.83/7.56 


10.28/12.02 



Fig. 19 shows the aircraft 10-g turn trajectory plotted in three 
dimensions for a 175 yards/second air speed. As compared to the 100 
yards/second air speed trajectory shown in Fig. 20, the distance 
traveled down range is much greater, and the turning radius was much 
larger. However, within the 5 second interval of flight, the air- 
craft in Fig. 19 did not climb as high as that in Fig. 20. This is 
because the third maneuver (turn up) did not start until approximately 
3 seconds versus 2 seconds for the 100 yards/second trajectory. 
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Figure 19 

Three Dimensional Aircraft Trajectory, 
Velocity = 175 yards/second 
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Figure 20 

Three Dimensional Trajectory, 
Velocity = 100 yards/second 
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In Fig. 21 the radial error and the measurement noise in yards 
is plotted versus time for the close-in maneuver. The radial error 
is equal to the true radial distance minus the estimated radial dis- 
tance in yards. Gaussian noise with zero mean and a standard deviation 
proportional to the square root of the estimated radial range variance 
was the measurement noise added to the true range. For this example 
the noise has a 1 sigma standard deviation of 2.0 yards. The very 
large peak error at 0.8 seconds is due to the aircraft flying very 
nearly "overhead". This then reverses the sign of the relative 
velocity vector. 

The velocity reversal is shown in Fig. 22 where X(2) is the Kalman 
estimator velocity state, and VR2D is the first difference of range 
divided by the sampling interval. The slow time constant of X(2) is 
due to the lack of a velocity measurement. A more accurate velocity 
estimate is shown in Fig. 23 for the distant turn maneuver. There is 
no velocity reversal in this trajectory. 

Fig. 24 shows the range acceleration state estimate, X(3), and 
the measured range acceleration, Y(3). The two accelerations track 
quite closely and describe two phenomenon: 

1. An initialization. 
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Figure 21 

Range Error and Measurement Noise 
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Figure 22 

Estimated Range Velocity and Differential Range Divided 
by the Sampling Interval, Close-In Maneuver 
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Figure 23 

Estimated Range Velicity and Differential Range 
Divided by the Sampling Interval, Distant Maneuver 
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Figure 24 

Range Acceleration Estimate and Range Acceleration Measurement 
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2. Cyclic acceleration matching the period of the third 
part of the maneuver. 

The acceleration state estimate in Fig. 25 with no acceleration 
measurement does not show this periodic acceleration, indicating the 
uncertainty present in the state estimate. 

This section applied the correlation attitude measurement tech- 
nique to the radar tracking problem. The attitude is related to 
acceleration and therefore a useful measurement in predicting target 
acceleration. Recursive trajectory equations were used in a computer 
simulation that estimated aircraft range. The simulation demonstrated 
that the Kalman Filter was more accurate in tracking a target with 
varying trajectories when the predicted acceleration was used as in 
addition to the noisy range measurements. 
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Figure 25 

Acceleration Estimate Without Range Acceleration Measurement 
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VII. CONCLUSIONS 



Manual optical image comparisons have been the primary method 
of measuring attitude of three-dimensional objects. Section V lists 
the general two-dimensional Fourier transform equations that are 
applied to an images edge data and a stored model of the solid ob- 
ject. The attitude of a three-dimensional object is estimated by 
using a three-axis quadratic interpolation to sequentially maximize 
the correlation output as a function of the attitude, position and 
size variables. Moreover, this measurement or estimate is equivalent 
to the maximum output of a bank of matched filters. 

This concept was then applied to realistic Kalman Filter track- 
ing simulations and shown to give improvement over conventional target 
state estimators that do not use the attitude measurement as a pre- 
dictor of target acceleration. 

Since the necessary attitude measurement is within the electro- 
optical tracker and computer state-of-the-art, the use of this esti- 
mator concept can improve present tracker system capabilities. 
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APPENDIX A 

TWO-DIMENSIONAL DISCRETE CORRELATION 

This appendix demonstrates that the two-dimensional discrete 
correlation transform (31), (32) is equivalent in both the time domain 
and the frequency domain. 

N-1 N-1 

c(k,l) = z z x(i,j) h(k + 1, 1 +j) (31) 

i=0 j=0 



C(m,n) = X*(m,n) H(m,n) (32) 

where * = complex conjugate 

Equations (31) and (32) are verified by performing the Fourier Transform 
and conjugation operations indicated. 
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The term in brackets in equation (A-2) is equal to; 



MN 



(M N) = 1 , if m - p and n = q 



0 , if m p and n q 

due to orthogonality of i and j. Therefore; 



2n 



c(k,l) = 



MN 



1 M-1 N-1 . ^ .• 

I I X*(m,n) H(m,n) e ^ 

m=0 n=0 



2n 



M 



n 1 



(A-3) 



= 2 ^ 



X*(m,n) H(m,n) 
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APPENDIX B 



COMPUTER PROGRAMS 

This Appendix is composed of the listings of the Fortran programs 
,nd subroutines for the attitude measurement and target tracking 
.imulations. The programs were run on a Univac 1110 computer with 
I Fortran compiler and a FLECS translator (Fortran Language with 
;xtended Control Structures), The graphical outputs were plotted 
ising DISSPLA programs that generated COM-80 crt images. 



I 
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C 1. 



1 / I c N b I fU w H ( 1 (, i C H ) . « M, ( 1 0 3 ! *. ) 

D I *1 fc u b I 0 N b ( 1 2 o ) . 1 1» V ( 1 c ' ) , { I ) 

Ol.ttNSIc-N AZAfo'.cLA(6^,KULLA<o'>,CP'A*f2”) 

0 I." fct. b I ON I (. t:0u ) 

COmklEX*o A(16ioA) 

COl'>l^U^ (•hh.WH.A.CflAX.Z 

DO in J = 1,2U 
C^A*(J)=-1.0c^0 

111 continue 

1 PL 0 r = 1 

I P « I f, T = 1 

1 hEAI/ (b*1u) AZ«EL«j\v>'LL*LAbT 
n F 0 A r A T ( ) 

aKI Tt(Oi1u1) AZ •£L*ROLL*LAST 
1J1 FOKE-rtT(1X,3Fy.2,lA) 

IF( LaST.EQ. 1 ) IPLUT=3 

L 

C READ ARRAY SIZE IN PO-ERS OF 2 

C READ NU^t;ER OF AHRAVS TO Bt*'. COhRELATEO fMAX OF 20> 

C 

R£AD(b.5)N<,JJ 
5 FORfAK^lS) 

READlb.S) NA 
N 1 = 2 * «K K 
n2=2**JJ 
n3= 1 

KK 2 =KK ♦ 1 

J J2 =J J + 1 

N12 =2»M 

n22=2*N2 

NH = M/2 

NJ“N2/2*»S 

N V = 1 * N 2 

NVA=n12*n22 

NC = (M*K2 )/ 2‘.M 

N ( 1 ) = KK 2 

?^^2 > = J J 2 
(« ( 3 ) = 0 
PI = 3. U159 
WRITE (6. 15) n1,n2 

15 FOHNAT(' IrjPUT ARRA'' onENilONb = <',Iif'X'.l3. '■>'./) 

C 

C READ IN TPE Data Array 

c 

WRITEfG«ii) 

DO 11 N=1 .NC 

K 1= ( K -1 ) » iA+ 1 

n2A=r«2‘. 

KEAC<5,7)nH(I),I=KlfKcA) 

••Rl lc(6*9)Al • (wr(I ) *I“N1 

11 continue 
^ F OH f AT ( 4 *, F i . L) 



66 



r\> 



FOk.-mtI' I ' .2^ I '-h ( I ) ') ./ / ) 

9 FOt<('AT(1X,l5,24(lX,F4.0)) 

C 

C OOUllE size of UaT« AKKAV 

1 . «><H = Input vtCTOR, whk = output vector, nv = input oincnsiov 

c 

call VD OUtJL ( «H , «iHR . nv ,NV <•) 

c 

C PLOT DATA ARkAY 

c 

c 

C TRANSFER THE Wh(K) DATA TO THE COMPLEX A MATRIX 

C 

00 i 12 = 1 .NZi 
DO 3 1 1 =1 ,N1 2 
IWH= (I2-1)*n 1*IV 
A(IWH) = WHR(It,H) 

3 continue 
C 

C COMPUTE THE DISCRETE FOURIER TRANSFORM OF THE luH(<) DATA 

C 

CALL HARM (A ,M , INV ,1> ,1 . IPERR ) 

C 

C READ IN THE SECOND DATA ARRAY 

C 

C 

C ROTATE MODEL IN AZf’-UTri 
C 

NL= 0 
NA1 =1 
OELAL =5. 

AZA ( 1 ) =AZ 
AZA (4)=AZ-DELAL 
AZA<3) = AZ-» DELAL 
290 CONTINUE 

DO 3C0 J=NA1 ,NA 
C 

C GENERATE IMAGE aOUNOAhX FROM 3 -0 I M E N T I ON A L MODEL 
C -HR = CmRD = OUTPUT VECTOR (INCYM X INCYM) 

C 

CALL SURFAC{AZA(J).EL.ROLL.YMAX.ZMAX,YMlN,iMIN,J.IPLOT.IPRlVT) 

C 

C CALCULATE CORRELATION BETWEEN INPUT IMAGE AND GENERATED IMAGE 30UnD» 
C 

G37 FOR mat ( / , 5X , ' J = ' , , 2 , 3 X , ' C M A x =',Fl7.(b,/) 

CALL C ALCOR (KK 2 . J J ,NV . AZ A ( J ) . E L . K OL L . J , N 1 2 , N2 2 ) 

30o CONTINUE 
NL= NL ♦ I 

IF(NL.GT,6) oO to 3202 
wRITc(6,3lO) 

31L FORMAT!/, 5X,'J',3X,'AZA(J)'.7X,'CMAX{J)', /I 

DO 311 J=1,3 

wRITcTe,312'' J,AZA(J),CMAX(J) 

311 CONTINUE 






it 




H'Ki-MT(sx.i2,3x.Fd.2,3*.f11.t) 

c 

L Ttsr FCn ^AXl^UN CORncLATIf'N. IF CI-AxCD < Ci*<AXf?) OH C^AX<j> 
C CoFREUATt FT ADblTIO.i^L FUIi^TS 
C 



I F ( L n A X ( 1 ) . L T . C .'1 A A i 


i ) ) 


GO 


T C 


3141 


I F ( C i>. A X n ' . L T . C M A X ( 


^ ^ 


hU 


TO 


31 A2 


GO TO 3144 










i 1 4 1 I F ( C F! A X ( 2 ) . L T . C rt A X i 


3 ) ) 


Gu 


T 0 


3 1 42 



AZA f 3) = A2 AM ) 

AZA(1)=A2A(2) 

AZA (<:) = A2 A(2) -DELAl 
O-A X ( 3 ) = C !»iA X ( 1 ) 

CMAXTDsCF'AXTS) 

C f" A X ( d ) = - 1 , E 1 0 

N A 1 = ^; 

N As L 

GO TO 2V0 
-• *1 H ^ CON TI NUE 

AZA (2) = AZA(1 ) 

AZA(1)=AZA<3) 

aza<3)=azaF3)*0ELAl 

CF'AX(2)=CFAX(1) 

C»^A X ( 1 ) =C."AX ( 3) 

CfA X ( 3) =-1 .E 10 

N A1 = j 
NA = 3 

bO TO 2VG 
i 1 A*, CON TINUE 

C CALCULATE MOST PRUoABLE AZMUTh FOR INPUT FhAME BY »^U0RATIC INTERF 

C 

AZST=(C«AX(2)-CiAAX(3))*0ELAL /(2.*CMAX<3)»Et*C>1AX(Z)-4.*CMAX(1)) 
IAZST=IFIX(AZST*.5)*IFIXTAZA(1)^ 

222c F0RF'AT(ZA6) 

PRIM 3 2C.AZST, lAZiT 

32o F OR F ATf / , 5X , ' AZ ST = ' , F 1 2 . 5 , 5x , ' I A Z S T s'. 14,/) 

C 

C REDUCE INTERVAL OF INTERPOLATION IF AdS. VALUl OF AZST > 1. 

C 

I F ( At S ( AZ ST ) .L T . 1 . u ) GO TO 320^ 

DEL AL = 2 tU 

NA 1 s2 
N As 3 

AZA(2)*AZA(1)-uELAl 
AZA f3)=AZAf1)*DELAL 
CMAX(2)=~1»E10 
CMAx( 3)=-1.E10 
GO TO 29u 
CONTINUE 
AZ= FLOAT! lAZ ST) 
c 

C RuTATE MODEL IN ElEVAIION 
C 



AP 



DEL tPS = !> . 

N L = 

c A X ( 7 > = C !•' A X ( 1 ) 

f, A3 =<i 

N A9 = V 

ELA (1 ) =tL 

ELA <^:) = EL-DELEl'J 

cLA(i)=£L*DELEPS 

CONTINUE 

DO AoO J = f.Aa#NAV 
jy= J-t 
L 

C GtNEPATc If"AbE oOUNOArtV FnOM 3 -D I i*i E N T I ON A L hOi/EL 
C WhR = CmhD = OUTPUT WtCTOR (INCTH X iNCY'i) 

C 

call SUhFAC(A2,ELA(J«),ROLL.YMAX.ZI*AX,Tf<'lN,ZWlN.J.lPLOT.IPRINT) 

C 

C CALCULATE CORRELATION bETwEEN INPUT IMAGE AND GENERATcD IMAoE 30UNDA 
C 

call CALCOR (KK2 .J J 2 .NV . A Z . EL A ( J H ) . R OL L . J . M 2 . N 2 ^ ) 

•.GG CONTINUE 
NL = NL ♦ 1 

IFTNL.GT.6) GO TO ‘.2C2 
»JRI Tt(6,A10) 

a1u F0RNAT(/,3x,'J',3a.'ELA(JI*)'.7X.'C«AX(J)'./) 

DO 411 J=7,9 
JM= J-c 

bRITE(0.312) J,tLA(JM).CMAX(J) 

411 CONTINUE 
C 

C TEST FOR MAXIMUM CORRELATION. IF CMAX(7) < CMAX(b) OR CMA*(9) 

C CORRELATE AT ADDITIONAL POINTS 
C 

I F < CMAX f 7) .L T .CMAX i 8> ' GO TO 4U1 
I F ( CMAX (7 ) .L T .CMAX ( 9) ) GO TO 4142 
GO TO 4144 

4141 I F ( CMAX (6 ) .L T .CMAX (9) ) GO TO 4142 
ELA (3) = ELA(1 ) 

ELA(1)=ELA(2) 

ELA(2)=£LA(2)-DcLEpS 

CMAXT9)=CMAX<7) 

CMAX (7) =CMAX (b) 

CMAX(b)=-1.ElO 
N Ab =0 
N A9 =0 
GO TO 390 
41 a2 CONTINUE 

ELA (2) = EL A(1 ) 

ELA < 1 ) = EL A(3) 

ELA (3)=ELA(3) + uELEPS 
C»"AX(b)=CMAX(7) 

CMA X (7) =CMAX (9 ) 

C’'AX(9)=-1.E10 
N A 9 = 9 
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I'V A :3 - 9 

b 0 T u i y G 
‘•lAs CCNTINUc 

c 

u CALCULATE >^OST PROcAdLt tLEVATION FOR INr^UT FRAVt c T liuORATIC I T 

ELST=(CFAX(5)-ct«AX<9))*C'cLEPS/<2,*Cf'Ax(9)*i,*CMA»fc1-4.*CHA*(1)) 
IELST = IFIX(ELST+.5>-»IFIX(ELA(1)) 

PRIr.T 321,cLST,IELi7 

12 ] FOHFAT(/,5x,'ELbT - ' , F 1 2 . 5 . Sx . M E LS T ='. 1 *,,/) 

C 

C RtDUCE INTERVAL OF INTERPOLATION IF At)S. vALUt OF ELST > 1. 

C 

I F ( Ao S ( EL ST ) ,L r , 1 . u ) 00 TO 4202 
0 EL EPS = 2 . 

NAd=d i 

NA? = 9 

ELA (2) = EL A(1 )-l>ELEPS 
ELA (i) = EL AM ) ♦OELtPS 
CHA X ( 8) = -1 , ElC 
C KA A ( 9 ) = -1 , £ 1 U 
00 TO 390 

continue 

EL=FLOAT(IELST) 

C 

C ROTATE MODEL IN ROLL 
C 

DEL RhO=5 . 

NL=L 

CMAX(13) = C?>1AX(7) 

ROL LA ( 1 ) = ROLL 
ROL LA( 2 )=ROLL-DELRhO 
ROLLA(3)=ROLL*DcLRrO 
NA1 A=U 
NA1 5=15 
49C CONTINUE 

DO 5C0 J=NA14,NA15 
J'-= J-12 
C 

C GENERATE IMAGE dOUNDAnY FhOH 3 -D I M E N T I ON A L hOuEL 
C WHR = CARD = OUTPUT VtCTOH (INCYM X INCYM) 

C 

call SURF«C(AZ.EL,R0LLm(JH),YMAX,2WAX.YMIN,ZMIN.J,IPL0T.IPkINT) 

C 

C CALCULATE CORRELATION BETWEEN INPUT IMAGE AND GENERATED ImAvjE 3 0UNDA 
C 

call CALCCR<KK2.Jj.l.NV ,AZ.EL.hOLLA(JM).J.Nl2.N22) 

50U CONTINUE 
NL = NL ♦ 1 

I F < NL .0 T . 6 ) GO TO i 202 
wRI T£ (6 .5 10) 

510 FORMAT(/,5x.'J',3x,'ROLLA(JK)'.7X.'CMAX(J) './) 

DO 511 J=15. 15 
J J -1 2 



rvj 



Wf'l Tt(c,31j) J.t\OLLA(JM).Ci«AX (J) 

i "! 1 continue 

C 

C TtST EOh I^AXI^UM COKrttLATION. IF CKAX(13) < Ci-AX(14) OH CC.A*(15) 

C C0HkEL><TE at AOdIFIONmL F-OIisTS 

L 

IF(C^’AX(1i).LT,Cf^A-i(lA)) GO TO 5141 
IF((.MAX(13 ).LT,CHiAa(15)) GO TO 51*.2 
GO TO :>144 

^141 IF(Cl*iAxn4).LT.CNiAAn5>' GO TO 514^ 

RCLLA(i)=KDLLA(1) 
hOL L 4 ( 1 ) = fvOLL A ( t ) 

ROLLAf2'=ROLLA(l)-uELf<hO 
CrMXl15)=C''. AX(13) 

CF-AX(131=C"Aa(1*,) 

CN-A X ( U ) = -1 . 1 1U i 
N A 1 4 = 1 4 
NA1 b=14 
GO Tu 490 
5 14c CON TINut 

KOL L A f 2 1 =ROL L A ( 1 > 

KOL LA ( 1 ) = KOLL A (3) 

ROLLA(3)=ROLLA(i) + LiELhHO 
C-AX(14) = C«iAX(13) 

C!«AXM3' = C''AAn3) 

Cf-AXCnis-Uclu 
NA1 4=1 5 
NA1 5=15 
GO TO 49C 
5144 CONTINUE 
C 

C CALCULATE MOST RRCoAciLE ROLL FOR INPUT FkAME 3y -UORaTIC INTE^P 

C 

ROLLST=(Cl'AX(l4)-t.^AX(15))*OELRHO/(2.*CF'AXU5) + 2.*CiMAX(lH) 

*-4. »CNAX( D) 

IfiOLSr=IFIX(RuLLST*.5)*IFIX(ROLLA( 1)) 

PRINT 323.ROLLST. IROLST 

3 FORKATC /. 5X, 'ROLLS r = ' , F 1 2 . 5 . 5 X . ' IROLST =',I4./) 

C 

C REDUCE INTERVAL OF iNiERPOLATlON IF AuS. XALUc OF ROLLST > 1. 

C 

I F ( AbS ( ROLLS T ) ,LT , 1 .0) GO TO 52''2 
DELRPO=2. 

NA1 4= 14 
NA1 5=15 

R0LLA(2)=R0LLA(1)-uELkH0 
RCLLA(3) = hOLLAC 1)*^ELR»*0 
CPA>Cl4)=-1.b1U 
CPAX(15)=-1.Elu 

GO To 4 V tj 

520c CONTINUE 
I PL OT = 2 

I F ( L A S T .E ) GO To 1 
STOP 



?1 



S L h f A c 



iUnf «C kOTaIES a STcKEO i-0 IMc.x T I on AL model In i AnIS A z r> u T n . t L E V A r - 
iOi\ Tht EhOJECTtU IE’AlE OF The I'^OOcL ON A ELANt SUkFaCc IS OuAnTIZcO 
ONTO M (32>3^) CKIO. OUTPUT ApkAY = 'CAku' 

S Lc ROUTINE SUkFAC(AZ.cL,ROLL .rf AXl.ZN'AXi .YMINl .ZMIN T.JA. IPLOT, 

, IPk INT ) 

PAmmi*'ETER NCS = 71«k 1 = 2*NCS 
p aAAMETER irJCYM =5i5 
? AHAMETE R I2E=C. I£=1 
pAmAMETcR NI =5 .nO = o. N0P = 7 

DiMtNSION AKfAX(rNCYM),P(3,-,,NCS).T(3,5.NCS).XN0RM(NCS).YK(2.NCS). 

A Y ^OK^ (NC S ) , Z NORM (NCS ) , AKM IN ( IN C YM ) . oKMAXdNCYM ), 

a bXMlNdNCYM ) 

D iMc NS ION CARO ( l63oA ) 

D I MEN SION YL(10),2L(10) 

C CMMON CARO 
C CMMON/ F AC 0 A T/ P 
Z ERU=U.O 

0 N£= 1 , 

L CTR=0 

D RC= ACOS Y-1 . ) / loC . 
l-ENERAL preperation of data 

AZR=AZ*LRC 
E LR = EL*D RC 
R OLL h = RuLL*0 RC 
S AZ=S1N(AZR) 

C AZ = CUS < AZR) 

S EL= SIN ( ELR> 

C EL=COS (ELR) 

S RO= S IN ( POLL R ) 

C hO= COS (ROLL R ) 

Y MA X = - 1 U . E 1 0 

Y MlN = 1(J. ElO 

Z RAX=-1U.E10 
ZMlN=10.Eir 

1 e22 F ORMATC 1X.I5 .12F 10.2) 
hOTA T E T AR6E T 

00 10 K=1,NCS 

0 0 10 J = 1 .A 

T (1,J,x;)=P(1 .j,k)*LEL*CAZ-P12.J.K)«CFL*SAZ*P(3.J .<) ‘SEL 
T <2.J .K)=P(1 .J ,KY*<CRu*sAZ*SRO«SEL‘CaZ )*p(2.J.<) *<C rO*CAZ-ShO*ScL* 

lSAZ)-P(ji.J.N)«SRO*CEL c-., 

T(i.J.X)=P(1.J.N)«(SR0*SAZ-uK0*SEL»CAZ)»P(c.J.K)*(Ct<0*SEL«S«Z 

1 ♦ S RO* C AZ ) ♦ P ( 2 . J . K ) * C RO *C t L 

1 F(T(2,J.K'.oT.y-'Ax' YMAX = T(2,J»K) 
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1 t?<Itx=T(2.J.K) 

I KT(i.J.»<).tT.Z^^AX) ZNaX = Ti;3.J.IO 

I KI<3,J,K’>.LT.2I'IN' Z'*IN=I<3.J.K) 

') c CNUUUt 

Ut FIRST YM — and vALUcS 

L CTR=LC TR* 1 
If(JA.Nh.l) OU ro 
Y f> A X 1 = Y t-t A X 
T H N 1 = Y n I N 
7 Kha 1 = Z F AX 
Z I* I N 1 = Z iM I N 
2 C CM I NUE 



,I>’>^UlE SURFACE NORMALS 



5 NORM<)=(T(Ll.X)-T(2.2.K))*(T(3,1,t^)-T(3.A.K))-(Tl2. 

, NO«n(K) = (T(3.1.Kl-T(3.2.IO).<T(1.1.lO-I<1.' ■'>>>-<’' • 

1 T ( 1 . A . X ) ) * < T ( 2 . 1 ) “T I i. . 2 . K ) ) 

15 CONTINUE 



OI^PU TE lYNAX AN 0 IZMAX 



L CTK=LC T R* 1 

R ZF‘AX=AbS(ZM AX-ZMN) 
R YfAAX = AfaSf XWAX-Yf-IN> 
' 1 ZKAX = RZF'AX 

I YhAX = RY>AAX 

l’ 

IrREPA re plot PAG E 



.1C 



I F ( 1 PLO T . E Q • C ) bC TO 19 
C ALL TMPPLT( 0 ) 

A XLL TH= AK.AX 1 ( f Z P A X - Z F» i N 
axloTh = aint(axlotfi/2») a 

X CNTR = AINT((YNAX + YPIN)/2) 

Y CNTR = AlNT((ZNAXtZPIN)/2) 



) , (YN AX 

10 



-Y f" I N 



X CKb = xCNTR-A XLOTH 

Y lRO = YCNTR-A XLbTP 
X cNO=xCixTR + AXLOTF 

Y tND=YCNTRAAXLGTH 

C ALL abNPL f - 1 ) 

C ALL NOd RO R 
C ALL PEIGPTC ,D 

E NCOCe ( e 0. 21 t . L Al-EL) AZ 

CALL T I T L E YL Ab EL . 1 ^ • 

C nLL GRm F ( XO I- b , 5 . X END . 
C K L L H A h K E R ( '- ) 

CALL SCLPIC(».jc) ^ 

f LhnAT < ' 



.EL.RULL 

YuRG»3«YE>xD) 



,F0 .1 . iX 



EL=' .F Y. 



X ) 



. 1 



R0LL= ' 



1 . < )- 
1.K) 

1 , X ) 



,F 6 . 1 



I w C Crjl If<L>L 
t 

ClIAhT IhLili FOR Z-KAiTtR 
t 

n i,C y i’ = f L 0 A T ( I N c Y rO 

0 LY=kY1'>aX/R><CYN 

D LZA = hZ|AAX/f(NcYh 

-.1 Fch."'AT(' PLT1. PL T > . ' i;6X . ' YMN . ' , 1GX , 'Z I'll n' . 1 • 'DL Y . ' , 1 cX . 

U OL^A • ^jXf IZP'AA* • 1 a* IYNAA •/) 

-t F o«>'AT ( a F 16 . 5, Z r. , / ) 

K J*AA =u 
K r Ih = L 

1 cc U = 1 • 

I iEcL = 1 

« hlTfc (6 , 1222 ) LCTrt 
1) C 1 CO I Z = 1 . lr,CYf 
Z tAK=IZ -OLZA+ZMIN 

MND Y INTEhSECTIONS 

L \.e= C 

A KKA X(IZ) = -1C'.E1G 
A Kf-IKCIZ ) = 10 .E1C 

2 51 FCRJ"AT(' Y, ZBARt A X f* I W ( I Z ) t A K M AX ( I Z ) t 0 Y y" I N t Z Y M I N ( IN C Y ) , D Y M A X , 
.Z YTcf«P.HAX(IKCY).',2QX,'IZ.lf IFLGtIAKFUj(lNCV).INCY'./) 

DC aG I = 1 ,NC S 
Y k(1 .1 ) = YMIN-100. 

I F (X NORh! ( I ) . LE .0 . ) 60 TO 4J 
N b K = 0 

00 25 X=1.A 
K E = K + 1 

1 F(k ,EQ .4) K E = 1 

I F<At3S < T f 3 ,K , I )- T< 3,K E , I ) ) .L E. 1 .E-6 ) 60 TO 25 

Y=(T(2.A.E.I)-T(2.K,I))*(ZeAK-T(3.K,l))/(T(3,Kt.l>-T<3,X,I))* 

1 T ( 2 .K . I ) 

X = (T(1.KE,I)-T(1,X.I))*(ZeAn-T(3.K,l))/(T(3.KE,l)-T (3,X.I))* 

A T( 1 , K , I ) 

I NC Y = ( Y - YM IN ) / OL Y 

15 intersection klTHlN cfOUNDS OF TFE SURFACE 

[ F(Y«LTiT(2fRfI)»AND»Y«LT»T(2*KE»II1 60 TO 25 
IF(Y.LT.T(2.K,1).AN0.Y.LT.T(2.KE.1)) 60 TO 25 

IF(ZcAR,LT.TCi.K.I).AND.ZdAR.LT.T(3.<E.I)) 60 TO 25 
I FfZbAR»6T.TY'i»x*i'»AN0*^BAr(.GT.T(3«KE.I') dO TO 25 
1 F (Y .LE . A«PA X ( I i. ) ) 60 TO 1 V :> 

A KF‘A X { I Z ) = Y 

IvS I F ( Y ,6E , AKMI N ( I Z > ) 60 Tu 19c 
A KF'iSf IZ) = Y 
1V6 C CNT INUE 
2C C CM INU £ 
i5 ,C CNT iNUc 
4 0 C C.N T I N U E 
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I 
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Ll NCIbt To AK>‘--{li) 



^LO I I-OI r. T S YL ( I ) 

IFUHlOT.EQ.v' UU to c 5 
^ L ( I ) = Z b A R 
Z L C) = Z w' A h 

»hEN(AKrAX(Ii;),LT.RLT1,OH.ANrA*(IZ).6T.HLTw) 

. ONLtbSFA<hiNTlZ'».LT.PLn.G(<.AM-'INfIZ).UiT.PLT2' 
. . Y L { 1 ) = /^^^ I N( 12 ) 

, . K I N = K (• I N ♦ 1 

. . CALL CUkVECYL.ZL. 1.-T) 

. ...fin 

...FIN ‘ 

E Lit 

. -H£N(An>fAIN(IZ).LT.PLT1.0n.AKMlN(IZ).(jT.PLT2) 

. . VL Y 1 T= AKI^AX ( IZ > 

. . K:R;AX=IC^AX + 1 

. . CALL C Uh VE ( YL . ZL . 1.-1) 

. ...FIN 

. ELSE 

. . YL ( 1) = AAMA X ( 12 ) 

. . Yl(2)=AKMIN(IZ) 

. . Ki" AX=< .XAX 4 1 

. . CALL CURVE (YL.ZL , 2.- 1 ) 

. ...FIN 

. . . F I Ti 
0*^ C ONT iNUt 
0 5 C Clx T iNUc 
iL'O C UNT INUt 

L LTR=LCTR4l 
J MN=0 
J t- A X = G 

stakt index for Y-RASTER 

L CAx =0 
L h I N = u 

wRITE(6,1222) LCTR 
0 C i CU I = 1 , INC YF' 

Y EAn = I « C'L Y ■* Y i«. I N 

FIND I INTERSECTIONS 

3nF'Ax(I ) = -1».E1w 
B M«I N ( I ) = 10 .E 1C 
OOaCu J=1.\Cb 

IF<XNOHF’YjY,Lt,(j.) 00 TO AO>j 
DO wLG ^•1.A 
K E = x 4 1 

IF(x.t0.4)<E = 1 

IF(«tb<TY2.K,j)-T(^i.KE.J))*LE.1.F-0) bO TO Z"0 



AKh- -( 1 ) = l-uK I 20M AL LIMtTS FOn SEi^U^nTIAl l SCAN 

ukm--(j)=vertical limits For ScCiUential y scan 



( CENTER Ii«AcE I 1 F- RESPECT TO FlhST IWAOE ANO 
I FUST IM SlAN image faOTTOM TO TOP (I) AND LEFT 

t DUCKET E Intervals. 

. CONVERT MA>, AND MN. EOoE COORDINATES TO OKID 
( 

D E Y K' I N = V N I N - t r- 1 N 1 
0 EYMAX=irF'AX-YhAXl 
0ELY = ( 0 E VMI N-*DE YM AX Y / 2 . 

Y MN=YMIN1-»0ELY 
D LYl=AcS(Y!^AXl-YMlNl)/RNCYM 
0 c2MlN=2MIN-ZrtIN 1 
DEZ«AX = 2F'AX-ZMAXl ‘ 

0EL2=( D EZMI N+ OE 2MAX ) / 2 . 

2 MIN = ZM1N1 + D ELZ 

0 LZAl = AcSTZ«Ax1-2MINl )/KNCYM 
I NCTh=0 

u 0 105 0 1 = 1 . INCYP! 

F X=(AKMAX(I)-ifMIN)/OLYl 

1 X=l FI X ( FX Y 

G X=(AaMIN(I)-YMIN)/0LY1 

1 G=I F I X (GX ) 

2 A character flag 

I K F = G 



CONTRAST PIXEL FLAG 
I tF= : 

00 IGoO J=1»INCYM 

F Y=(bXMAX(J)-2MN)/DLZA1 
J F=I F I X ( F Y ) 

G Y=(bKMIN(J) -ZMIN)/DLZA1 
J G=I F I X ( GY ) 

IR =Y1-1)*INCYM*J 

1 F(I R.Nc .0) GO TO 103 j 
I R = 24 

I RF= 1 

1 u30 C ONT INUE 

IF(IX.Nc.J) GO TO 1033 
CaH D ( I H ) = 0N E 
I EF= 1 

GO TO 1 044 
U33 C OnT INUE 

CA«D(IR) = ZERO 
IF(IG.NE.J) GO TO Ij3g 
CAkO ( IR Y = ONE 
I EF= 1 

I NC T R = I N C T R ♦ 1 
GO [0 1 LA4 
1 i.it C CNT INUE 

CARO ( I R ) = ZE RU 
I F(J F .NE .1 ) Go TO 1 j4 J 
CAR C ( I R ) =0N E 
1 EF= ^ 



USE 6 R I ki 
TO RIGHT 

integers 



IfJTERVALS 
( J ) AT 



0 F 
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A T( i , ^ . J ) 

t iS ItiTchScCTIOK wIlhiN aOUNUi Ot- Tht SUKFACt 

If(ZF»LT,T(i«KfJ),A'i)D»Zf»LT«r(j'. «k.c»J)) oO TO 30 j 
lMZF.GJ.Tf3,K,J^.AN0.Zf.GT.T(3,»E,J') bO TO 3 00 
IMYLA‘<.GT.T(i:.K.J),AND.yBArt.GT.T(Z.KE.J)) bO TO 30 u 
lF(Yt;rtA.LT.T(4.K..J).AND.YBAh.LT.T(2.Kc.J)) oO TO 30 o 
IFfiF.Lt.bKMAXfl )) GO TO 
= Ki- A X ( I ) = 2F 
L i*~A=LAiHX+1 
C t.-,T iNUc 

IFf^f.Gt.BKnKJM )) GO TO c 
) = ZF 

L I-IN = Li>' 1 1 

c55 C L.NT INUE 
3Cr C ONT INUt 
A L 0 C uN T I h U E 

ADO KANuCM NOISE TO 3KM--(Ii) 



FOOT F-OInTS 

IFMPLOT.EG.C) GO TO ASS 

Y L(Z ) = Yd AR 

Y L(1 ) = Yti Ah 

W PEN (tjKci AX ( I ) ,L T .PL T3 .OR .fa^^'A X ( I ).6T.PLT4) 

. U NL E S S < BK M N f I ) . L T . P L T 3 . OR . B K h I N T I ).GT.PLTA’» 

. . ZL(1) = rhMN(I ) 

. . JMN = JMN+1 

. . CALL CURVE (YL.ZL. 1 , - 1 ) 

. . . . F I r< 

. ..FIN 
E Lib 

. i.Y'cN (BKhlNCl ).LT»PLT3.UR.t;XiMlN(I ).bT.PLT>.) 
. . ZL f 1 > = EKF AX ( I ) 

. . JKAX=JFAX-» 1 

. . CALL CUF.VETYL.ZL. 1.-T) 

. . . . F I fj 

. c L S £ 

. . 2 L ( 1 ) = t h ^ A X ( I ) 

. . 2l < 2) = LKhIN( I ) 

. . J i'l A X = J F A X ♦ 1 

. . J I N = J t- I N •* 1 

. . CmLL CUKVE(YL.ZL. b.-i) 

. ...FIN 

. ..Fix 

~c‘ C ur.T if^Uc 
- 0 ‘ C C F. f I 111 U c 
CLNfluUt 
L LTF=LCrh+1 
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a 





I 



M 

i 









UC fu 

1 •« n c c.N r I Mj t 

CArvD^ Ik^ = ZE hO 
1 t- ( J 6 . N t . I ) b 0 TO 1 '1 <, i 
C An 0 ( I K ) = ON £ 

I tf = 1 

0 0 TO 1 C44 
114 3 C b-NT INUt 

CARD(Ifv)=ZEKU 
1 < h4 C CNT InUE 
1 vCG C CNT INUC 
1 (.5 0 C CNT INUt 

L CTh=LC T 1 

1 fdPRiNT.Ea .c) tiO TO 

P hlr^T 2oCfAZ t&LfPuLL 

C nUL END PL (0 ) i 

L CTk=LCTR+ 1 
MC C ON T INUE 

L CTh=LCTR* 1 
R ETURn 

iuC FORhATdPl. 5X,'AZIC*UTH A NS L E = ' . f 8 . 2 , / . 5 X . 

A 'ELEVATION AN(iLE = ',F8.<l./ , JX, 'ROLL A N b L E = ' , F 8 . 2 . / / J 
c NO 
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iC'3'^ul.T INE VuC'ucLfARlN.i^hOT.lNLi'^.NUO) 

C 

(. VtiCUhJL CONVEkTS Ai^ iNhUT (NXN) AhHAY ItjTO A (^NA2 n) AhRAy fcJT rtDDING 
V. L'ACkgrluno data, 

C AkIN=INPbT ARRmY, MkOT=CUTPUT AdUhY, INOrl=itvPUT OlKtMlON 

C 

UIMlNSION a R I N ( I NUi'’- ) . APOT(NOO) 

NOD = a*1Mj(^ 

FO^‘ = fLOAT(I^u^') 

RDM = bUft T ( FDM ) 

1DR = IFI*(KD^') 

IDN£' = i:*ItiK 

Zf FOHF AKZX ,'IUR 15,/) 

L 

c jsKOw Index for Output array i=colunn index 
c 

DO ^uU J=1.IDR2 

I F ( J ,GT .1 DR ) GO TO 27C 

DC £00 1=1, IDR 

IJ= IDR2*( J-1 )*I 

IJlC = I0R*(J-1)-»I 

AP0T(IJ)=ARIN(1JK) 

£..1 continue 

c 

C FILL remainder uF R0»* WITH ;CER0S 
C 

DO 250 1 = 1, IDR 
IJ= IdR2*( J-1 )-HOR^ I 
A R 0 T ( I J ) = 0 . 

25 u CONTINUE 
GO To 300 
27U CONTINUE 
C 

C FILL remainder of R0«Z WITP ZEROS 
C 

00 2oO I=1.IDR2 
IJ= IdR2*< J-1 ) ♦ 1 
A R 0 T ( IJ ) = 0 . 

(^CLj CONTINUE 
iOu CONTINUE 
RETURN 
END 
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C THIS SUBROUTINE CALCULATES THE CORRELATION BETWEEN 

C the a VECTORTOFT COEFFICIENTS) AND THE B VECTURfOFT COEFFICIENTS) 
C 

SUBROUTINE C AL C OR ( KK 2 , J J 2 , NV , A Z , EL . R OL L . J .N 1 2 . N 2 2 ) 

DIMENSION WH(16384),WHR(1638A) 

DIMENSION SM28),INV(123),Mf3) 

DIMENSION CMAX(20) 

DIMENSION MAXX ( 20) .MAX Y (20 ) 

DIMENSION W( 12b. 12&) 

C0MPLEX*8 Arl6i<34) ,e(l63b4) 

COMMON bHR.WH.A.CMAX 
M(1 ) = KK2 

M(2 )=JJ 2 » 

M(3 )=0 

2037 F0RNAT(/,5X,'J = ' . I 2 . 3 X , 'C M AX =',F17.6./) 

63 FORMAT( 5X ,'n12 ='.I4,'N22 =',I4.'NW =',I4) 

C 

C DOUBLE SIZE OF DATA ARRAY 

C HHR = INPUT VECTOR. WH = OUTPUT VECTOR. NV = INPUT DIMENSION 

C 

NV4 =4*N V 

CALL VDOUBL (WHR ,WH ,NV ,NV4) 

C 

C PLOT DATA ARRAY 

C 

C TRANSFER THE WH <K ) DATA TO THE COMPLEX B MATRIX 

C 

NDA TAsQ 

DO 104 12=1. n22 
DO 103 11=1. N12 
IWH= (I2-1)*N12+I1 
B(IWH)=WH (IWH) 

IF< WH(IWH) .6T.0.0) NOATA = ND AT A* 1 
103 CONTINUE 

IWHM=lWH-9 

WRITE(6.72) (WH ( I ) ♦ I=IWHM. IWH ) 

104 CONTINUE 

URITE(6.208) NDATA 
20o FORMAT!' NDATA ='.17) 

I F( NOATA.EQ .0) STOP 
C 

C COMPUTE THE DISCRETE FOURIER TRANSFORM OF THE wH (K) DATA 

C 

71 FOR MAT( 16 ( 1 X , F7 .2) ) 

72 FORMAT (10(1X.E12.5)) 

73 F0RMAT(7(2X.I5) ) 

call harm (8 . M , I NV . S . 1 . 1 PERR ) 



C 

C 



COMPLEX multiply 



ijv. 




I 

I 




e(I 1)=A(I1)*C0NJ6(b(I1)) 

IF(n .LT.N22) WRITE(O.IOI) I 1 , A ( II ) , 8 ( 1 1 ) 

220 CONTINUE 
C 

C inverse FOURIER 

C 

CALL HARM (8 ,M , INV . S ,-1 . IPERR) 

C 

C CALCULATE THE MAGNITUDE OF THE FOURIER COMPONENTS 

C AND STORE ONLY THE POS FREQ VALUES OF ii(I1,l2) 

C 

DO 2U I 2=1, N22 
DO 3U 11=1. N12 
IWH= (I 2-1)*Nl2*I1 
M n 1 , 12 ) = CAB S (8 ( IWH ) ) 

I F( W< II . I 2) .LT.CMAX (J ) ) 60 TO 30 
CMAX<J)=W(I1,I2) . 

MAX X ( J ) =I 1 
HAXY(J)=I2 
30 CONTINUE 
C 

C PRINT THE OFT COEFFICIENTS 

C 

IF( MAXX (J ) .GT.3) 60 TO 40 
MM3 = 1 
MM2 =2 
MM1=3 
MAX X ( J ) =4 
MP1 =5 
MP2=6 
MP3 = 7 
MP4=tt 
GO TO 42 
40 CONTINUE 
N12M=N1 2-4 

I F ( MAXX ( J ) .6 T .N 12M ) GO TO 41 
MM3 =MAX X ( J ) -3 
MM2 =MAX X ( J ) -2 
HM1 =MAX X < J ) -1 
MP1 =MAXX ( J) + 1 
MP2=MAXX( J)+2 
MP3=MAX X ( J ) ♦ 3 
MP4=MAXX( J)*4 
GO TO 42 
41 CONTINUE 
MM3=N12-7 
MM2=Nl2-6 
MM1 =M2-5 
MAX X ( J ) =Nl2-4 
MP1 =n12-3 
MP2=n12-2 
MP3=N12-1 
MP4=M2 
42 CONTINUE 



TRANSFORM 
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hfil Tt(6.1t0) Mni,HW2,hM1,f'AXX rJ),»">^1,HP2,KK'3,f'P4 
1 :u FORr’AT(' K'.5X.'A('.I3,'.K)',7(6X.'A(' ,I3,',K)')) 

101 fORf-AT(2X.I4,8(lX,kl2.5)) 

00 61 12=1, N22 

WRITE(6.1U1)I2,(W(I1,I2).I1=MM3.MPA) 

6l CONTINUE 
RETURN 
END 
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c 

C KAL is a radar range ESTlf»ATOR USING A KAL^'AN FILTER. 

C THERE ARE 2 RANGE MEASUREMENTS AND A RADIAL ACCELERATION ESTIMATION 
C MEASUREMENT. 

C 

C MODEL 

C X(K*1) = H(0*x(lO+6(0*U(t() 

C Y(K+1) = S(K+1 )*X<K41 )>V(K^1 ) 

C KALMAN estimate 
C X'(K+1) = H(k)*X-HAT(K) 

C X-HAT(K*1) = * ' (K )-K (K^1 )* (S (K + 1 )*X ' (K )-T (x+1 ) ) 

C k(K + 1) = P ■' ( K ♦ 1 ) *S T (K ♦ 1 ) * ( S ( K ♦ 1 ) *P ' ( K^ 1 ) * S T ( K ■* 1 ) *R ( K 1 ) ) * • -1 
C P'(K + 1) = H(K)*P(K)*HT(K)-»6(K)*Q(K)*6T(K) 

C P(K + 1) = P'(K*1 )-K (K + 1 ) *S (K + 1 ) •P'Ck^ 1 ) 

C P(0) = G(0)*0(0)*GT(0) 

C T = TRANSPOSE 
C subroutines USED: INTCON 

C MANVR 

C FtKALI 

dimension X(4),H(4,4),6<4),S<3,4) ,P<4.4),T1(4,4),T2(4,4),T3(3), 
*TX(4),Y(3),R(3,3),E6<4.4) 

dimension Tpnoo) ,EP1 (100) ,XP1 (100) .RNOlSnOO) ,VP1 (100) ,VP2(100) , 
*RNOIS2(100),XTP(100).2TP(100) 
dimension API (100).AP2(100) .AP3(100) 

DIMENSION WHEN(2) ,TME0D(2) 

DIMENSION PAUV(3)*ANUV(3)«ANL(3),ANLR(3)tTPA(3) 
dimension 10PT(5).STAT(5> 
dimension IPAKdSO) 

COMMON XT,YT,ZT,VXT,VYT,VZT,VEL0C,T1AN6,T2AN6,DELR.IG,TRAT1, 
.DELT,NPRF,TRAT2,II,NPH,R90,I62.ICTS,P1.IR,IH,I0S 
COMMON/INTC/ A Z . E L . ROL I ♦ R 0L2 . R 0 L4 , A Z D E 6 , E L D EG , ROL 1 06 .R0L2D6. 
♦ROL406 ,T6S1 ,T6S2.RT.R0OTT»R20OTT, D T , T R A D 1 * TR A 0 ? . R A N 6T , 
•RANGTM,XR,YR,ZR.VXTM,VYTM,vZTM,VR20,AXT,AYT,AZT,AXT6,AYTG,A2T6, 

* ANA ,ANE »PAUV 
DATA IN,N,IT/3*4/ 

DATA IS,M1/2*3/ 
data IL.L/2*1/ 

DATA ISEED/42573/ 

DATA PI /3 .i 41 592654/ 

C 1/(S**2) model state VECTOR 

CALL FR80I0(' 6220 R.H. WILLIAMS 6262 ') 

Nl = 5 
N0 = 6 

TGS1=10. 

TGS2=10. 

0T=0.05 

NP = 31 

NPM=NP-1 

NPH=NP/10 

RES=0.5 

R90=PI/2. 

R270=3.*PI/2. 

6YO=10.7252 
NRUN = 0 
10 CONTINUE 



c 

C 3-0 TARGET TRAJECTORY DATA 
C ORIGIN AT INITIAL TARGET POSITION 
C RT = TRUE RANGE RELATIVE TO RADAR 
C ROOTT - TRUE VELOCITY RELATIVE TO RADAR 
C R?DOTT = TRUE ACCELERATION RELATIVE TO RADAR 
C VELOC = TARGET VELOCITY RELATIVE TO ORIGIN 
C trad- = TURN RADIUS OF T6S- 6. TURN 

c TRAT- = Turn rate of tgs- g, turn (pad/sec) 

C -R = RADAR location AT T(0) 

C -T = TRUE X«y,2 target POSITION 

C V-T =: TRUE X,Y,Z VELOCITY 

C A-T = TRUE X,Y,Z ACCELERATION 

C A-TG = TRUE X«V,Z ACCELERATION <G'S) 

C ANA = VERTICAL (A7) ACCELERATION RELATIVE TO TARGET ORIENTATION 
C ANE = HORIZONTAL (EL) ACCELERATION RELATIVE TO TARGET ORIENTATION 
C PAUV(I) = PITCH AXIS UNIT VECTOR 
C PAUV(1)=X 
C PAUVC2) = Z 

C PAUV(3) = Y 

C 

C ATC = TINE CONSTANT OF ACCELERATION 

C TCK = TRANSITION NATRIX ACCELERATION TIME CONSTANT 
ATC = 2. 

TC = ATC/OT 

TCK = ( 1 .-E XP (-TC ) ) 

C a = COVARIANCE MATRIX OF U (LXL) - INPUT 
C = ECU*UT] = (RMS RAN6E/SEC*SEC*SEC)**2 
R£AO(S,150) Q»RCV,ILAST 
00 115 10S=1,2 
IF(I0S.EQ,2) ot=o.i 
DO 118 IR=1,3 
00 116 IH=1,2 

c H = state transition matrix (NXN) 

DATA H/16*0./ 

C AGMTK = gamma acceleration TRANSITION MATRIX CONSTANT 
C ARTK = NORMAL RADIAL LOAD ACCEL. TRANS. MATRIX CONSTANT 
IF(IH.EQ.I) AGMTK=1.0 
IF(IH.E0.2) A6MTK=0.0 
IF(IH.Ea.3> AGMTK=0.5 
IF(IH.EQ.I) ARTK=0.0 
IF(1H.EQ.2) ARTX=1.0 
IF(IH.E0.3) ARTk=0.5 
H(1 ,1 ) = 1 . 

H(1 ,2)=0T 

H(1,3)=AGMTK*0T*DT/2. 

H(1 ,4)=ARTK*DT*DT/2. 

H(2,2) = 1 . 
h(2,3)=DT*AG"'TK 
h(2,4)=ARTK.0T 
h(3,3)=1. 

C ANLRK = RADIAL NORMAL LOAD ACCELERATION TRANSITION CONSTANT 
A NL RK =1 . 

H ( 4 ,4 ) =ANLRK 

c G = state noise input matrix (NXl) 
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300 F0RWAT(' Q =',F12.1,' R =',F9,3»' RT =',F10.1.' ROOTT = ' , 
*4(F7.1),I3) 

CALL OATE(WHEn) 

CALL TOFOAY (T»tEOD ) 

WRITE(N0.112) WHEN(1),WHEN(2),TMEOD(1)tTWEOO(2) 

112 format (2A6 ,2X ,2A6) 

WRITE (6,90) 

90 F0RMAT(' K',4X,'RT ',7x,'Y',3X. 

*4X, 'X (1)', 7X,'X(2)‘',2X, 'ERROR =RT-X(1)',2X, 'NOISE', 
•2X,'AXT6',2x,'AYT6',2x,'A2T6',2X,'2N0 OIFF-VELOC',/) 

WRITE(6,130) X ,RT , Y (1 ) ,X (1 ) , X (2) , EP1 ( 1),RN0ISd) 

C 

C MAIN LOOP 
C 

DO 100 11 = 1, NPW 

IP=II+1 

IPS*0 

IF(II.Ea.l) IPS=1 
PCTNS=0.7 ‘ 

RN0IS(II)*SQRT<R(1 ,1))*66N0F(ISEED)*PCTNS 

RNOIS2(II)=SQRT(R(2,2))*GSNOF(ISEED)*PCTNS 

OELR=VELOC*OElT*FlOAT(NPRF) 

C MANVR CALCULATES TRUE POSITIONS (XT,YT,2T) AND TRUE VELOCITIES (VXT, 
C VYT,V2T) TARGET STARTS S ORIGIN 3 T'O, AND FLIES TOWARD RADAR FIXED 
C LOCATION RT ON X AXIS 

C X = DOWN RANGE, REFERENCED TO TARGET S T =0 

C Y = UP RANGE, REFERENCED TO TARGET 3 T =0 

C 2 = OFF RANGE, REFERENCED TO TARGET 3 T =0 

CALL MANVR 

C V = total true VELOCITY RELATIVE TO COORD SYSTEM , T(0) 
V=SQRT(VXT*VXT+VYT=VYT+V2T*VZT) 

AXT=(VXT-VXTM)/(DELT* FLOAT (NPRF)) 

AYT=(VYT-VYTM)/(DELT* FLOAT (NPRF)) 
A2T=(V2T-VZTM)/(0ELT*FLOAT(NPRF)) 

C AT = total TRUE ACCELERATION RELATIVE TO COORD SYSTEM , T(0) 
AT=SQRT(AXT*AXT*AYT*AYT^A2T*A2T) 

AXTG=AXT/10,7252 

AYT6=AYT/10,7252 

A2TG=AZT/10.7252 

ANE=((AYT*VXT-VYT*AXT)*VXT-(A2T*VYT-V2T*AYT)*VZT)/(V*VH) 

ANA=(A2T*VZT-VZT*AXT)/VH 

C AH = horizontal ACCELERATION - RELATIVE TO COORD. SYS., T(0) 

AH = SQRT(AXT*AXT-*-A2T*AZT) 

ELE=ANE+COS (EL) 

ELN=SaRT(ANA*ANA+ELE*ELE) 

C VH = HORIZONTAL VELOCITY - RELATIVE TO COORD. SYS., T(0) 

VH=SQRT(VXT*VXT+V2T*V2T) 

A2=ATAN2(V2T,VXT) 

A2DEG=AZ*180./F>I 
EL=ATAN?(VYT , VH) 

El0£G=EL*180./PI 

anael=ana/ele 

C ICTS = coordinated turn switch 

IF(ICTS.EQ.I) ROL4=ATAN2(ANA,10,7252) 

R0L5 =ATAN2 (ANA,10.7252) 
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6(1 )=DT*OT*DT/A, 

G(2)=0T*DT/2. 

G(3)=0T 
6(4)= G (3 ) 

C R = covariance wATRIX of V (MlXMl) - t^EASUREHENT 
DATA R/9*0*/ 

R (1 ,1) = RCV 
R(2,2)= R(1 ,1) 

R(3,3)=100‘R(1,1) 

C S = OBSERVER MATRIX (mIXN) Y = $X 
DATA S/12*0./ 

S(1 ,1 )*1. 

S(2,1)=1. 

S(3,4)=1. 

C INTCON INITILIZES TRAJECTORY VARIABLES 
C 

CALL INTCON 

C TX = ITERATED EXTENDED THRESHOLDS 
TX ( 1 )*1 . 

TX (2)=TX (1 ) 

TX (3)*TX (2) 

TX (4)=TX(3) 

T = 0, 

LL*0 

INCP=0 

RN0IS(1)>SQRT(R(1t1>)*66N0F(ISEED) 

RN0IS2( 1)=SQRT(R(2,2))*66N0F(ISEED) 

C Y = INPUT OBSERVATION VECTOR (HiXl) 

C multiple MEASUREMENTS 
V(1 ) = RT4RNOIS (1 ) 

Y(2)=RT*RN0IS2(1) 

Y(3)*0. 

C PRESET estimator INITIAL STATE 
C X(1) > RANGE(K) 

C X(2) = RANGE RATE(X) 

C X(3) = RANGE ACCELERATION 

C X(4) = ESTIMATED RADIAL NORMAL LOAD ACCELERATION 

XlOOFF= 1. 

X20k=1 .1 

X(1)=Y(1)^XlOOFF 
X (2)=X20X*RD0TT 

X (3)=0. 

X (4)=0. 

EP1 ( 1 )= RT-xn ) 

XP1(1)=X(1) 

VP1 (1 )=X(2) 

VP2 (1 )=RDOTT 
AP1 (1 )=X(3) 

AP2 (1 )=0. 

AP3 (1 )=0. 

FPMax=0, 

R0RX=RT 

NRUN=NRUN+1 

WR1TE(6»300) Q,RCV.RT,R00TT,A6MTX,ARTK*S(3.3)fNRUN 
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144 FOR«AT(' COOO. TURN S U ='I3,' ROL4 =',F9.3,' COORD. TURN ROLL 
, CV.3) 

ROLOE6=ROl4*180./PI 

ANAD0T=<(-VXT*AXT-VZT*A2T)*ANA)/(VH*VH) 

azoot=ana/vm 

eldot=ane /V 

R OLOOT = 1 0.7252* AN A00T/(1 0.7252*10. 7252 ♦ANA *ANA) 

C V-TM = DELATED TRUE X.Y.Z VELOCITY 

VXT«=VXT 
VYTM=VYT 
VZTM=VZT 

C -3 = X.Y.Z DIFFERENCE BET*,EEN RADAR AND TRUE TARC-ET POSITION 

X3=XR-XT 
Y3=YR-YT 
Z 3 = ZR-ZT 

C RAN6T = TRUE RANGE RANGTM = DELAYED TRUE RANGE 

RAN6TM=RANGT 

RAN6T=SQRT <X 3*X3-f.Y3*Y3*Z3*Z 3) 

VN0IS=RN01S (II) 

RANGN=RANGT*VN01S 

RANGN2=RANGT*RN0IS2(II) 

RE*(DR=AHOD(RANGNtRES) 

RE»»DR2 = AH0D(RAN6N2»RES) 

RES2=RES/2. 

C RANGE = MEASURED RANGE WITH NOISE AND QUANTIZED TO RESOLUTION 

RAN6E=R ANGN-REMDR 
RAN6E2=RANGN2-REMDR2 
RT=RAN6T 

C AWNL = MAGNITUDE OF NORMAL LOAD ACCELERATION ( G'S ) 

AL*1 3 . 

BE«-3 . 

GAM=0.5 

EPS=66N0F (ISFED) 

AMNL*AL*BE*EXP(6AM*EPS) 

C AZ = AZDEG 
C EL = -ELDE6 
C ROLL = -R0LDE6 
C PAUV(I) » X 
C PAUV(2) = Z 
C PAUV(3) = Y 

SAZ=SIN( AZ) 

CAZ»COS( AZ ) 

SEL=SIN (-EL) 

CEL»COS(- EL ) 

SRO = S IN (-R0L4 ) 

CR0=C0S(-R0L4) 

TPA(1)=PAUV(1)*CEL*CAZ-PAUV(2)*CEL*SAZ*PAUV(3)*SEL 

TPA(2)=PAUV(1)*(CRO*SAZ*SRO*SEL*CAZ)*PAUV(2)*(CRO*CA2-SRO*SEL*SA 

.-PAUV(3)*SR0*CEL 

TPA(3)=PAUV(1)*(SRO*SAZ-CRO*S£L*CAZ)*PAUV(2)*(CRO*SEL*SAZ*SRO*CA 

. *PAUV (3)*CR0*CEL 

C ANUV(I) = normal acceleration UNIT VECTOR 
C = (TPA(I) X V-T)/V (VECTOR CROSS PRODUCT) 

C (X ,Z ,Y) 

ANUV ( 1 ) =( (TPA (2)*V YT) - (TPA<3)*VZT))/V 
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ANUV (2 ) = ( (TPA (3)*VXT) - (TPA(1)«VYT ))/V 
ANUW(3)=( (TPA(1)*vZT) - (TPA(2)*VXT))/V 
C ANL(l) = NOPWAL load ACCELERATION 
ANL(1) = ANUV(1)*ANNL 
ANL(2) = A NUV ( ? ) * ANNL 
ANL(3) = ANUV(3)*A«NL 

C ANLP(l) = RADIAL normal LOAD A C C E L E R AT I ON ( 6'S ) 

C = PROJECTION (DOT SCALAR PRODUCT) OF ANLd) ON THE RADIAL 

C VECTOR FROM TARGET TO RADAR/RADIAL RANGE 

ANLR(1)=((ANL(1)*x3)+(ANL(2)*23)^(ANL(3)*y3))/RAnGT 
Y ( 1 ) GRANGE 
Y (?)=RANGE? 

Y ( 3 ) =ANLR (1 ) *GYD 

C VflPO = 1ST DIFFERANCE, RANGE VELOCITY VR?DD = DELAYED VR?D 
VR2D0=VR2D 

VR?D= (RANGT-R ANGTM) /(0ELT*FL0AT(NPRF)) 

C AR = RANGE ACCELERATION, 2ND DIFFERANCE 
ARM=AR 

AR=<VR2D“VR20D)/(DELT*FLOAT(NPRF)) 

' WRITE(6,133) RANGT,RAN6TN,VR2DD,VR20,AR,Y(3).X(3) , 

*X3,Y3,Z3 

C TEST RADIAL 2ND DIFFERENCE: AR, & LIMIT TO MAX, EXPECTED 
C RATE OF CHANGE: SQRT(Q) 

ADI F=ABS ( AR-ARM) 
eLlM=SQRT(Q)*DT 

IF(ADIF.GT.aLIM) AR=SIGN(QLIM,AR)+ARM 
C TEST FOR RE-INITILIZ ATION OF KALMAN FILTER 
SI6TH»2.0 

S1G1 = SIGTH*SQRT(R(1,D) 

SI62=SI6TH*SQRT(R(2,2)) 

SI63 = AMAX1 (SIGI ,SI62) 

PHAX*0. 

EEST=((Y(1)>Y(2))/2.)-XC1) 

IF(£EST .LT.S163> GO TO 44 

INCP=INCP^1 

IPS*1 

DO 42 JP*1,N 

PMAX=AMAX1(P(JP,JP),PMAX) 

42 CONTINUE 
PPCT = 1 .0 

DO 43 JP=1,N 

P( JP, JP)=PMAX*PPCT 

43 CONTINUE 

44 CONTINUE 

WRITE <6, 133) AZ,EL.R0L4,TPA(1),TPA(2) ,TPA(3),ANUV(1) ,ANUV(2), 
,ANUV(3),AMNL,ANLR(1)fAR 

CALL FTkAlI (K , X ,H ,6 , Y , S ,0 ,R , P , I N , I S , I L , N , M 1 , L , T I , T 2 , I T , T 3 , 
•TX,F6,IER,IPS) 

T*T*DT 

EP1 ( IP) = RT-X ( 1 ) 

XP1(IP)=X(1) 

VP1 <IP)=X(2) 

VP2(IP)=VR20 
API (IP)=X(3) 

AP2 (IP)=Afl 
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AP3(lP)=y (3) 

TP ( 1P)=T 
XTP(IP)=*T 
ZTP (1P)=ZT 

VPCT=X(2)*OT*100./EP1(IP) 

APCT=X(3)*DT*OT*50./EP1(lP) 

PCT» = T2(1 ,1 )*100./EP1 (IP) 

WRITE(6,134) 1I,VPCT,APCT,PCTH,T2(1»1),PMAX 
134 FORWAT (IX ,13,' PCT VELOV, ACC, CEAS = ' , 1 0 ( F 1 2 . 3 ) ) 
1F(ABS(EP1(1P)),6T,ABS(EPMAX)) eP«AX=£Pl(lP) 

I f ( ABS(EP1 ( IP) ) ,GT ,40. ) LL=LL+1 

1 F (F6 (1 ,1 ) ,GE .1 . ) GO TO 51 
WNS=FG(2,1)/((1.-FG(1,1))*DT) 
lF(teNS.LT,0.) GO TO 51 
WN=SQRT (WNS ) 

WNHZ=>iN/2,*P1 

0R=(F6(1,1)-DT'*FQ(2,1))/(2*(1.-FG(1,1))*0T*WN) 

194 F0RMAT(' NATURAL FREQ. (HZ) =',F9.2,' DA'^PING RATIO ='.F8.3) 

51 CONTINUE 

WRITE(6,130) K,RT,V(1),X(1),“X(2),EPl(IP),RN0IS(II),AXT6,AyTG,AZTi 
*VR2D 

100 continue 

C 

I0PT(1)=1 
lOPT (5)=1 
NPP=NP+1 

CALL BEIU6R(EP1,NPP,I0PT,STAT,IER) 

HRITE(6,102) STAT(1),STAT(5),EPMAX 
102 FORMAT! ' ARITHMETIC MEAN OF RANGE ERR =',F7.3, 

VARIANCE (AVE. RMS ERROR) OF RANGE ERROR =',F9.3, 

MAXIMUM RANGE ERR =',F9.3,/) 

WR1TE(6,104) LL.INCP 

104 FORMAT!/,' NO. OF RETURNS OUT OF GATE =',I4, 

•' NO. times P INCR. =',I4) 

I RPL0T=1 

IF(IRPLOT.EQ.O) GO TO 113 
CALL INTAXS 
CALL TAXANG(0.) 

CALL SIMPLX 

CALL TITLE!' KALMAN F I LTE R S' , 1 00 , 

* 'ERROR IN YARDSJ', 

. 100, 'time in SEC. s', 100, 6. ,8.) 

CALL XTICKS(5) 

CALL YT1CKS(5) 

X0R6=AMIN1(ARMIN(EP1,NP),ARMIN(RN0IS,NP)) 

XMAX6=AMAX1(ARMAX(EP1,NP),ARMAX(RN0IS,NP)) 

CALL GRAF (XORG, 'SCALE', XMAXG, 

* ARMIN(TP ,NP) , 'SCALE', AR'«AX(TP ,NP)) 

CALL CURVE(EP1, TP,NP,10) 

CALL CURVE (RNOI S ,TP ,NP , 10) 

CALL MESSAG('COV Q =$',100,1., b.9) 

CALL REALN0( Q,0, 'ABUT', 'ABUT') 

CALL MESSAG(' COV R *$',100, ' AB U T ' , ' A B U T ' ) 

call RFALN0(R (1 ,1 ) ,2,'ABUT' ,'A8uT') 
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CALL «ESSAG(' » =*■'♦ 100, ' APUT ' , 'AbUT ' ) 

CALL lNTNO(N(?UN,'AflUT','A8UT') 

CALL OATE(WHEN) 

CALL TOFDAy (TMFOD) 

CALL MESSA6(WMEN,10,1.,9.3) 

CALL MESSAGCTWEOOtlO.'ABUT'.'ABUT') 
lDU«Mt=LlNEST(lPAK,150,l30) 

CALL L1NES('ERR0R»',1PAK,1) 

CALL L1NES('N01SES',IPAk,2) 

CALL LEGENDCIPAK .2, 4. 0,2.0) 

CALL ENDPL(O) 

113 CONTINUE 
I VPL0T=0 

1 F (IR .EQ.3) 1VPL0T=1 
I F (IVPLOT.FQ.O) 60 TO 114 
CALL INTAXS 
CALL YAXAN6(0.) 

CALL SIMPL* 

CALL TITLEC' KALMAN F IlTE R S' , 1 00 , 

* ' VELOC IN YAROS/SECt', 

. 100, 'TIME IN SEC. $',100, 6, ,8.) 

CALL XTICKSC5) 

CALL YTICkS(5) 

X0R6=AMIN1(ARMIN(VP1,NP),ARMIN( VP2,NP)) 
XMAX6=AMAX1(ARMAX(VP1,NP),ARMAX( VP2,NPj) 

CALL GRAF (X0R6, 'SCALE', XMAX6, 

* arminCtp ,np ), 'scale ' ,ARMAX (TP ,NP)) 
CALL CURVECVPI, TP,NP,10) 

CALL CURVE(VP2,TP,NP,10) 

CALL MESSAGC'COV Q =$', 100,1., 8,9) 

CALL REALNOC 0 , 0 , ' A8 UT ' , ' ABUT ' ) 

CALL message' COV R =$',100, ' AB UT ' , ' AbUT ' ) 

CALL REALN0(R(1,1),2,'A8UT','ABUT') 

CALL MESSAGE' # =$', 100, ' ABUT ', 'ABUT ' ) 

CALL INTNOENRUN, 'abut', 'ABUT') 

CALL DATEEWHEN) 

CALL TOFOAY ETMEOO) 

CALL ME SSA6 EWHEN , 10 ,1 . , 9 . 3 ) 

CALL MESSA6ETME0D, 10, 'ABUT', 'ABUT') 
IOUMMt=LlNESTElPAK,150,b0) 

CALL L I NES E 'X E 2 ) S ' , IPAK , 1 ) 

CALL LINESE'VR2D$', IPAK ,2) 

CALL LEGENDEIPAK ,2, 5. 0,1,0) 

CALL ENOPLEO) 

114 CONTINUE 

I aplot=i 

1 F EIAPLOT.EQ.O) go to 1148 
CALL INTAXS 
CALL YAXANGEO.) 

CALL simply 

CALL TITLEE' KALMAN F 1 LT E R $ ' , 1 CO , 

« ' ACCEL IN YOS/SEC 'SECI' , 

. 100, 'time in SEC, $',100, 6. ,8.) 

CALL XTICKSE5) 

CALL YTICkSES) 



X0<?6 = AMIN1(ARM1N(AP1,NP) ,A»MIN( AP^.NP)) 

XMAX6 = A*«Ax1(AR»»AX(AP1,NP),aR>«aX( AP?,NP)) 

CALL 6R A F ( X 0R6 , 's C AL E ' , X M AX G » 

• ARi»I1N(TP ,NP ) , 'sc ale ' , a RI* ax (TP ,NP)) 
CALL CURVE(AP1, TP.NP.10) 

CALL CURVE (AP2.TP.NP.1C) 

CALL CURVE( AP3.TP.NP. 10) 

call MESSA6('C0V Q =»', 100,1.. 8.9) 

CALL REALNOC Q.O.'APUT'.'ABUT') 

CALL MESSAGC' COV R =$'.100, ' AB UT ' , 'ABUT ' ) 

CALL REALN0(R(1.1),2,'ABUT','ABUT') 

CALL MfSSAGC' » =S', 100, ' ABUT ' , 'AEUT ' ) 

CALL INTNO < NRUN , 'ABUT ', 'ABUT ' ) 

CALL DATE(WhEN) 

CALL TOFDAV (Tr»EOD ) 

CALL ME S S AG ( WHEN , 1 0 , 1 . , 9 . 3 ) 

CALL MESSA6(TME0D, 10, 'ABUT', 'ABUT') 

IDUMi«v = lINEST(1PAK,150,80) 
call LINES('x(3)S',1PAK,1) 

CALL LINES('AR $',IPAK,2) 

CALL LINES('Y(3)$',1PAK,3) 

CALL LEGENDdPAK ,3,A.0,?.0) 

CALL ENOPL(O) 

11A8 CONTINUE 

CALL INTAXS 
CALL YAXAN6(0.) 

CALL SIMPlX 

CALL TITLE ( 'TARGET GROUND PLOT S', 100, 

*' OFF RANGE (2) - YDS. S', 100, 

*' DOWN RANGE (X) - Y D S . S ' , 1 00 , 6 . , 8 . ) 

CALL XTICKS(5) 

CALL YTICKS(5) 

X0R6= AR«IN(ZTP,NP) 

XMAX6= ARMAX (ZTP ,NP) 

YMAXG=ARMAX(XTP,NP) 

CALL GRAF (X0R6, 'scale'. XMAX6, 

• ARMINCxTP ,NP) , 'scale' ,YMAX6) 

NPLS= 5 

CALL CURVECZTP.XTP.NP.NPLS) 

C DRAW ♦ AT RADAR LOCATION, XT=RT , 2T=0. 

CALL MARKER(3) 

CALL SCLP1C(2.) 

RDR Y=0. 

1 F ( RDR X .LE . YMAX6) CALL C U R V £ ( R D R Y , R D R X . 1 , 1 ) 

CALL RESET< 'SCLPI C') 

CALL RESET('MARKER') 

CALL MESSAGC'COV Q =S', 100,1., 8.9) 

CALL REALNO( 0 , 0 , ' A B U T ' , ' A faU T ' ) 

CALL MESSAGC' COV R =S',100, ' AB U T ' , ' A BUT ' ) 

CALL REALN0(R(1,1).2,'A8UT','ABuT') 

CALL MESSAGC' « =S', 100, ' AB UT ' . ' A B U T' ) 

CALL INTNOCNRUN, 'ABUT', 'ABUT') 

CALL OATE(WHEN) 

CALL TOFDAV ( THEOO ) 

CALL MESSa6(WHEN,10.1 . ,9. 3) 



^ • 



I 




CALL «ESSA6(TME00,10,'AbUT','ABUT') 

CALL ENOPL(O) 

116 continue 

118 CONTINUE 
115 continue 

I f < ILAST.NE .1 ) GO TO 10 
CALL OONEPL 

120 F0RI«AT(lX,IA,1X,f6.1,3(F6.2),6(F9,3),6(F7.3)> 
130 FOPWAT (IX, 13. 4 ( F10.1 ) ,7( F 10.3) ) 

133 FORf'AT (1 X , 1 4 F9 .2) 

150 FORI^AT () 

END 
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SUPROUTINE INTCON 

C INTCON INITILIZES TRAJECTORY VARIABLES 
DIMENSION PAUV(3) 

COWMON *T,YT.ZT»VXT.VVT»VZT,VEL0C,T1AN6,T?ANG,DELR,I6,TRAT1, 
.DELT,NPRF,TRAT?,II»NPH,R90,1G2,ICTS,PI,IR,IH,I0S 
COMMON/INTC/ AZ.ELtROLl ,R0L2,R0L4,AZDEG,ELDE6.R0L106 .ROL2D6, 
*R0LAD6 tTGS 1 ,T6S2 ♦RT.fiOOTT ,R2 dOTT , D T » TR A 0 1 ♦ T R A D 2 , R A NGT * 
*RAN6TM,XR,YR,ZRtVXTM»VYT«,VZTM»WR2DfAXT,AVT,AZT,AxT6,AYTG,AZT6, 

*ANA , ANE ,PAUV 
ITCS=1 
AZ=0. 

EL = 0. 

ROL1=0. 

ROL2=0. 

R0L4=0. 

AZDEG=0. 

EL0E6=0 . 

ROL106=0. 

R0L2DG=0. 

T1AN6=0. ( 

16*0 

162=0 

DELT=OT 

NPRF=1 

C 3-0 TARGET TRAJECTORY DATA 
C ORIGIN AT INITIAL TARGET POSITION 
C RT = TRUE RANGE RELATIVE TO RADAR 
C ROOTT = TRUE VELOCITY RELATIVE TO RADAR 
C R200TT = TRUE ACCELERATION RELATIVE TO RADAR 
IF(lR.Ea.l) RT=11275. 

IF(IR,E0.2) RT=1000. 

IF(IR.Ea.3) RT=150. 

IF(IOS.EQ.I) RDOTT=-175. 

IF(10S.EQ.2) ROOTT=-100. 

R200TT=0. 

C VELOC = TARGET VELOCITY RELATIVE TO ORIGIN 
VELOC=-ROOTT 

C TRAD- = TURN RADIUS OF T6S- G. TURN 

TRAD1*VELOC«VELOC/T6S1*10.7252 
TRA02=VELOC ‘VELOC /TSS2* 10. 7252 
C TRAT- = TURN RATE OF T6S- G. TURN (RAO/SEC) 

TRAT1=TGS1‘10.7252/VELOC 
TR A T2=T6S2*1 0.7252/ VELOC 
RANGTsRT 

RAN6TN=RAN6T-(RD0TT*0ELT*FL0AT(NPRF)) 

C -R = RADAR location AT T<0) 

XR=RT 

YR=0. 

ZR=0. 

C -T = TRUE X,V,Z TARGET POSITION 

XT=0. 

YT = 0. 

7T = 0. 

C V-T = TRUE X.YtZ VELOCITY 

VXT=VELOC 
VYT=0. 

VZT*0. 

VXTN=VXT 
VYTr=VYT 
VZ TW = V2 T 



VR?D=RD0TT 

C A-T = TRUE X,Y,2 acceleration 

A * T =0 . 

A YT=0. 

AZT=0. 

C A-TG = TRUE X,Y,Z ACCELERATION (6'S) 

A X T G • 0 • 

A YT6 = 0. 

AZTG=0. 

C ANA = VERTICAL (AZ) ACCELERATION RELATIVE TO TARGET ORIENTATION 
C ANE = HORIZONTAL (EL) ACCELERATION RELATIVE TO TARGET ORIENTATION 
ANA^O. 

ANE=0. 

C PAUV(I) = PITCH AXIS UNIT VECTOR 
C PAUVd) = X 
C PAUV(?) = 2 
C PAUV(3) = Y 

PAUV(1)=0. 

PAUV (2) =1 . 

PAUV(3)=0. 

END 
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FLECS ,S 



MVR ,f< .t^VKP , h'(/RT 



manvr calculates true positions IXT,YT,ZT) anp true velocities (VXT, 

WYT.VZT) TARGET STARTS il ORIGIN S T’C, ANU FLIES TOWARD RADAR FIXED d 
LOCATION RT ON X AXIS 

X : DOWN RANGEf REFERENCED TO TARGET 3 T -C 

Y - uP RANGE, referenced TO TARGET a T rr 

Z r OFF RANGE, REFERENCED TO TARGET a T rO 

SUBHOUTINE PANVR 

CCMVON XT,YT,ZT,VXT,VYT,VZT,VEL0C,T1AN6,T2ANG,0£LR,IG,TRAT:, 
.0£LT,NPRF,TRAT?,II,NPh,KRr,lG2,ICTS,Pl,IC',IH,I0S 
COPMON/INTC/ AZ,EL,RCLl,ROL2,RGLRtAZDEG,FLD£G,PCLlDG ,RCL20G, 

♦ R0L4DG , TGSl ,T6S? ,RT ,RDOTT ,R200T T , D T , TR A D 1 , T R AQ2 , » ANG T , 
♦RANGTH,XK,YR,ZR,VXTM,VYTm,v2TP,VP2P*AXT,AYT,AZT,axTG,AYTG,AZT6, 

«AN A , ANE tPAUV 

»P0L4 ALONG X-AXIS ROTATES PITCH VECTOR CLOCKWISE. (PROP X,Y,Z - 
C,D,1 TO G,-A,P ) 

RP r ROLL RATE OF AIRCRAFT (DEG/SEC) 

RPR : POLL RATE OF AIRCRAFT (PAD/SEC) 
iCTS = COORDINATED TURN SWITCH 
DTrDELT 
R R ~ c E • 

PRRrPfi *PI/ 18C, 

CONDITIONAL 
. (II.LT.NPH) 

. . INITIAL TRAJECTORY 

. . ITCSrO 

. . X TzxT ■» DELR 

. . 'measured* aircraft POLL, PRECEEDING MANUVFR 

. . R0L‘* = R0L4 ♦ ( RRRidOT ) 

. . . .FIN 

. (T1AN6.LT.R90) 

. . SECOND TRAJECTORY 

. . R0L4=P0L4 ♦ ( RPR*DT ) 

. . IF (R0L4.GT. J .55) 

. . . R0L421.E5 

. . ...FIN 

. . IF(IG.CT.ift) 

. . . ICTS=0 

. . . 'ME/SUREO' aircraft roll, PPEGFECING MANUVEP 

, . . ROL4=ROL4-(RRR*DT ) 

. . ...FIN 

. . IG=IG>1 

. . TIANG : A2 ANGLE OF TARGET 

, , TiaNG:TRaTi<‘D£LT*FLOAT( NPRF)«I6 

. . XT =X T ♦CELR* COS ( T 1 ANC ) 

. . ZT=ZT ♦DELR*S IN ( T IANG ) 

. . GR9n VXT : L, VZT : VELOC 

. . VXT=VEL0C>>C0S ( T I ANG) 

. . VZT rvEL0C<‘S IN ( T 1 f-NG ) 

. ...fin 

. ( .TRUE . ) 
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. . iHlt-C TKSjE’CTG^Y 

. . •('■tiSUfLr’ A I fi CRAFT POLL 

. . IF ( K 0 L 4 . t T . C . C ) 

. . . R0L4:PCL4-(PRR«DT ) 

• • 

. . IF «ROL 4 ,L T .. .b ) 

. . . KOL**rr.o 

F IN' 

. . Ib2rIG2+l 

. . Ti: At.G : TP A T?<=OLLT»FLOA T ( NPRF ) *IG 

. . YT=YT4r.'ELR*SlN« T2ANG) 

. . ZT:2T4DlLP*C0S< T2ANG) 

. . VYT::vFLOC*S IN t T2ANG) 

. . VZ T=VELOC«COS ( T: A N6 ) 

• • • • F I \ 

• • • F I N 
RETL-ON 
END 
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.CS VERSION P 



C SUPROuTJNE FTk*LI 

c-f Tkah s 

c * 

f 

c Function 

C USA6F 

c 

C PARAMETERS K 

C 

C 

C 

C 

C X 

C 
C 
C 

C H 

C 

C G 

C Y 

C 

C S 

C Q 

C 

C R 

C 

C P 

C 

c 

c 

C IN 

C 

C IS 

c 

C IL 

C 

C N 

C 

C Ml 

c 

C L 

C 

C T1 

C 

C T? 

C 

C IT 

C 

C T3 

C TX 

C lER 

C 

C 

C 

C 



(XtX,H,6,Y,S,Q,R,P,IN,IS,IL,N,«1,L,T1,T2, 



IT,T3,TX , lER) 

- ITERATED XALMAN FILTERING 

- CALL FTk A IE ( X , X ,H , 6 , Y , S , 0 ,R ,P , I N , I S , I L . N ,M1 ,l 

T1 ,T2 ♦ IT,T3 ,TX, lER) 

* INPUT STEP COUNTER* WHEN X~0? 

VECTOR X SHOULD CONTAIN THE PRIOR ESTIMATE 
OF THE MEAN OF X, AND THE PROGRAM 
CALCULATES THE ESTIMATED VARIANCE Of 
X AS P=6Q6-TR ANSPOSE AT STEP 0. 

- INPUT/OUTPUT VECTOR OF LENGTH N. ON INPUT, 

X IS THE state vector AT STEP K, AND ON 
OUTPUT, X CONTAINS THE ESTIMATED STATE 
VECTOR AT STEP X*1. 

- INPUT matrix of DIMENSION N BY N. H IS THE 

transition matrix at STEP K, 

- INPUT MATRIX OF DIMENSION N BY L AT STEP K. 

- INPUT OBSERVATION VECTOR OF LENGTH Ml AT 

STEP K+1, 

- INPUT matrix of dimension Ml BY N AT STEP 

- INPUT COVARIANCE MATRIX OF U, OF DIMENSION L 

AT STEP K, 

- INPUT COVARIANCE MATRIX OF V, OF DIMENSION 

Ml BY Ml at STEP 

- INPUT/OUTPUT MATRIX Of DIMENSION N BY N. 

ON INPUT, P IS the variance MATRIX OF X 
AT step X. ON OUTPUT, P IS THE ESTIMATED 
VARIANCE MATRIX Of X AT STEP X-M. 

- FIRST DIMENSION OF H,6* AND P AS DIMENSIONED 

IN calling program. 

- FIRST DIMENSION OF S AND R AS DIMENSIONED IN 

CALLING PROGRAM. 

- FIRST DIMENSION OF 0 AS DIMENSIONED IN 

CALLING PROGRAM. 

- INPUT, SEE DESCRIPTIONS OF X,H,G,M,P. 

N MUST BE greater THAN 0. 

- INPUT, SEE DESCRIPTIONS OF Y,S,R,T3. 

Ml must be GREATER THAN 0, 

- INPUT. SEE DESCRIPTIONS OF G,Q. 

L MUST BE greater Than 0, 

- WORK MATRIX OF dimension NM BY nM, WHERE 

NM IS THE MAXIMUM OF N AND Ml, 

- WORK MATRIX OF DIMENSION N* BY NML, WHERE 

NML IS THE MAXIMUM OF NM AND L. 

- FIRST DIMENSION OF T1 AND T2 AS 

dimensioned in calling program. 

- WORK VECTOR OF LENGTH Ml. 

- INPUT, ITERATED EXTENDED THRESHOLD 

- ERROR parameter. I£R = 0 IMPLIES NO ERROR. 
TERMINAL ERROR = 1?e>N 

N = 1 INDICATES ONF OF IN, IS, IL, OR IT 
IS TOO small, OR that ONE OF N, *»1, 

OR L IS NOT A POSITIVE INTEGER. 






c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



c 



c 

c 



c 



c 

c 



c 

c 



^ = ? indicates I'ATf.I* inversion failed, 
PRECISION - SINGLE 

REG'O IR'SL routines - L £ 0 T 1 F , LUO A T F , L UE L '' F , Uf fi T S T , V » UL F F , V UL F P 
LANGUAGE - FORTRAN 



LATEST REVISION - HAY 10. 1V79 

hoDEL 

X(K 41 ) = H (K ) * X (< ) + 6 (K ) *u( K) 

Y(X,^1) = S ( K •» 1 ) • X ( 1 ) ♦ V ( K ♦ 1 ) 



SUBROUTINE FTxALI 



(K,X,H,6,T,S.Q,R,P,IN,IS,1L. 

N,M1,L,T1,T2,IT,T:*,Tx,F6,IEP,IPS) 



DIMENSION X(1),H(1N,1),G(IN,1),¥(1),S<IS,1),(S(IL.1). 

• R ( IS ,1 ) ,P( 1N,1) ,71( IT ,1 ) ,T2( IT , 1 ) ,T3 ( 1 ) 

dimension Xx(20), XP( 20) . E X (20) , Tx < in) , T6 (20) »HNL (2C ) , 

*EXM(20) ,FG(IT,1) 
dimension 75(20), T6(10, 10) 
data 10/0/ , 11 /I /, 120/20/ 

IF (IN.GE.N .AND. IS.GE.M1 .AND. IL.GE.L 

* .AND. (IT.GE.N .or. IT.6E.M1) .AND. N.GT.O 
^ .AND. Ml.GT.O .AND. L.GT.O) GO TO 5 

lER = 129 
60 TO 9000 
5 lER = 0 



calculate P if K = ZERO 

P(0) = 6(0)*0(0)*6T(0) 

IF (X .NE. 0) 60 TO 10 

CALL VMULFF ( 6, Q, N, L, L , I N , I L , T 2 , I T , I E R ) 

CALL VMULFP (72, G, N, L, N.IT.IN, P.IN.IER) 
IF(IPS.EQ.O) 60 TO 166 

WRITE(6,160) <P(1,J),J=1,IN),(TX(I),I=1,IN) 

DO 165 1=2, IN 

WR1TE(6,161) (P( I ,J ) , J = 1 , IN) 

165 CONTINUE 

166 CONTINUE 

160 FORMAT(' P(0) = ' , 1 0 ( E 1 0 . A ) ) 

161 FORM AT (S X , 1 U ( E 1 0. A ) ) 

10 CONTINUE 

SAVE INITIAL X 
DO 12 1=1, N 



XX ( I )=X ( 1 ) 



12 CONTINUE 



CALCULATE X-PRIME AT STEP X •» 1 



X'(X+1) = H(X)*X-HAT(X) 

CALL VMULFF (H, X, N, N , I 1 , I N , 1 N , T1 , 1 T , 1 E R ) 

DO 15 I = 1 , N 

XP(I) = T1(I,1) 

X ( I ) =XP ( 1 ) 

15 CONTINUE 

1 F ( I PS . £0 . 1 ) hRITE(6.17G) ( X p ( 1 ) , I = 1 , N ) 

170 FOBMAK' X-PRIME = H.X = ' , 1 0 ( F 1 0 . 2 ) ) 

CALCULATE P-PRIHE AT STEP K*^ 



p'(Xfl) = H(X ) *P (X ) *HT (X )♦& (X ) *0 (X ) *6T (X ) 

CALL VMULFF ( H, P, N, N, N , 1 N . 1 N , T 1 , I T , I £ R ) 



CALL VWULFP (T1, H, N, N, N,IT,1N, P,IN,1£R) 

CALL VULFF ( 6, Q, N, L, L , I N , I L ♦ T 2 , I T , I E R ) 

CALL Vf*ULFP (T2, 6, L, N , I T , I N , T 1 , I T , I £ R ) 

00 25 I = 1 ,N 
DO 20 J = 1,N 

P(IfJ) = P(I,J) ♦ T1(I,J) 

20 CONTINUE 
25 continue 

I F ( IPS.EQ .0) 60 TO 186 
WRITE(6,180) (P(1,J),J=1,IN) 

00 185 1 = 2 , IN 

WRITE(6,1S1) (P(1,J),J=1,IN) 

1S5 CONTINUE 

186 continue 

180 FORHATC' P-PRIME = ' , 1 0 < E 1 0 . 4 ) ) 

161 FORMAT(11X,10(F10.4)) 

c calculate matrix k at step tc^l 

C K(K + 1) = P'<K4l)*ST(K+1)*(S<K+1)AP'(K*1)*ST(K*1)^R<t«: + 1))**-1 
CALL VNULFP ( S, P,Ml, N, N , I S , 1 N , T 2 , I T , I £ R ) 

IF(IPS.EQ.O) 60 TO 28 
00 27 1=1, Ml 

WRITE(6,ie7) I,(T2(I,JJ),JJ=1,M1) 

27 CONTINUE 

28 CONTINUE 

CALL VMULFP (T2, S,M1, N , Ml , I T , I S , T 1 , I T , I £ R ) 

00 35 I = 1 ,M1 

DO 30 J = 1,M1 

T1(I,J> = T1(I,J) ♦ R(J,I) 

30 CONTINUE 

IF(IPS.EQ.I) nRIT£(6,187) I , ( T 1 ( I , J J ) , J J ’ 1 ,M 1 ) 

35 continue 

187 FORMATC' S.P.ST^R(',I1 ,',J) = T1 = ' , 1 0 ( F 1 0 , 2 ) ) 

CALL LEQT1FCT1, N ,M 1 , 1 T , T 2 , I 0 ,T 3 , 1 E R ) 

IF (lER ,EQ. 0) 60 TO 40 
lER = 130 
60 TO 9000 
40 00 50 I = 1,M1 

DO 45 J = 1 ,N 

T1(J,I) = T2(1,J) 

F6(J ,I )=T1 (J ,1) 

45 CONTINUE 
50 CONTINUE 

IFdPS.EO.O) 6a TO 194 
WRITE(6,190) (T1 ( 1 , I ) , 1 = 1 ,M1 ) 

DO 193 1=2, IN 

WRITE(6,191) (T1(I,J),J=1,M1) 

193 CONTINUE 

194 CONTINUE 

190 FORMATE' K = T1 = ' , 1 0 < F 1 0 . 2 ) ) 

191 FORMAT(10X,10(F10,2)) 

c CALCULATE X-HAT AT STEP K+1 

C XALMAN ESTIMATE 

C X-HAT(K + 1) = X'<K)-K(K^1)<*(S(K-M)*X'(X)-T(X41)) 

KL = 0 

DO 53 I = 1 , IN 
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f X ( 1 ) = 0 . 

F X -X ( I ) = 1 F 9 

53 continue 
5A CONTINUE 

00 542 1 = 1, IN 

T4(1)=XP(1)-x(1) 

542 CONTINUE 

C OBS. FRR. = S (X'-X)-Y4hNL = T5 

55 CALL V»!ULFf ( S, X,M1, N , 1 1 , 1 S , 1 N , T 3 , 1 S , 1 E « ) 

CALL VWULff ( S,T4,M1, N , 1 1 , I S » 1 N , T 5 , 1 S , I E R ) 

CALL VHULFFC S, X,M1, N , 1 1 . 1 S , I N , HN L , I ? C , 1 E R ) 

00 60 I = 1,M1 

T3(I) = T3(I) - Y(l) 

60 CONTINUE 

IFdPS.EQ.1) «RITE(6,195) ( T 3 ( 1 ) , I = 1 , P 1 ) , ( X ( I ) , 1 = 1 , N ) 
195 FORWAK" OBS ERR = S.X-Y =',12(F9.2)> 

C T1 = GAIN K 

CALL V«ULFF (T1.T3, N , W 1 , 1 1 , I T , I S , T 2 , I T , I E R ) 

00 65 I = 1,N 

X(I) = XP(1) - T2(l,1> 

T6(I ,1)=T2(I .1) 

65 CONTINUE 

IF(IPS.EQ.I) WR1TE(6.200) ( X ( I ) . I = 1 , N ) . ( T 2 ( 1 , 1 ) . I = 1 , N ) 
200 FORMATC' X-HAT = ' , 1 0 ( F 1 0 . 2 ) ) 

69 CONTINUE 

c calculate p at step X -»1 

C P(X-»1) = P' (X ♦! ) -K (K ♦ 1 ) *S (X + 1 ) *tp' (K+ 1 ) 

CALL VNULFF (T1, S. N,N1, N . I T , I S ♦ T 2 . i T , I E R ) 

CALL VNULFF (T2, P, N, N, N , I T , I N ,T 1 , I T , I E R ) 

DO 75 I = 1,N 
T2 (I ,1)=T6(I ,1) 

00 70 J * 1,N 

P(I,J) = P(1,J) - T1(1,J) 

70 continue 
75 CONTINUE 

IF(IPS.EQ.O) 60 TO 2l6 
WRITE(6,210) (P(1 ,J) ,J = 1 ,IN) 

00 21 5 1 = 2 , IN 

\xRITE(6,211 ) (P( 1 , J ) , J = 1 , IN) 

215 CONTINUE 

216 CONTINUE 
60 TO 9005 

9000 CONTINUE 

CALL UERTST ( 1 E fi . 6H F T K A LM ) 

9005 RETURN 

120 FORX'AT(1x,14,1X,F6.1,10(F11.2)) 

130 F0RNAT(1X,I5,12(F10.3)) 

140 F0Ri«AT(2X , 1 3 ,1 2 (F10.2) ) 

210 FORKATC' P =',10(F10.2)) 

211 FORNATC 5X , 10( F 10 .2) ) 

END 



100 



BIBLIOGRAPHY 



1. By: Infared Countermeasures Branch, Infared Signatures: Theory, 

Techniques, and Equipment . NWC Tech. Publ. 4902, 1971. 

2. S. A. Dudani, An Experimental Study of Moment Methods for Automatic 
Identification of Three Dimensional Objects from TV Images , Ohio St. 
Univ., PhD. Diss., 1973. 

3. R. L. Pio, "Euler Angle Transforms", IEEE Trans, on Auto. Control, 
Vol. AC-11, no. 4, pp. 707-715, Oct. 1966. 

4. E. H. Ehling, Range Instrumentation , Prentice-Hall Inc., Englewood 
Cliffs, N. J., 1967. 

5. A. V. Oppenheim and R. W. Schafer, Digital Signal Processing , 
Prentice-Hall, Inc., Englewood Cliffs, N. J., 1975. 

6. W. C. Effinger, Fortran Subroutine to Calculate Projected Surface 
Areas of Rotated Three Dimensional Objects, NWC Memo Reg. 4055-34-73, 
1973. 

7. A. Rosenfeld, "A Non-Linear Edge Detection Technique", Proc. of the 
IEEE , pp. 814-816, May 1970. 

8. R. L. Fox, Optimization Methods for Engineering Design , Addison- 
Wesley, Reading, Mass., 1971. 

9. E. 0. Brigham, The Fast Fourier Transform , Prentice-Hall, Englewood 
Cliffs, N. J., 1974. 

10. J. D. Kendrick, and others, "Estimation of Aircraft Target Motion 
Using Pattern Recognition Orientation Measurements," presented at 
the IEEE Decision and Control Conference, San Diego, CA., December 
1978. 

11. M. D. Levine, "Feature Extraction: A Survey", Proceedings of the 

IEEE, Vol. 57, pp. 1391-1407, August 1969. 



12. I. E. Abdou and W. K. Pratt, "Quantative Design and Evaluation of 
Enhancement/Thresholding Edge Detectors", Proceedings of the IEEE , 
Vol. 67, pp. 753-763, May 1979. 

13. R. L. T. Hampton, A Guided Tour of Kalman Filtering , NWC Technical 
Note 4070-58-75, 1975. 

14. By: International Mathematical and Statistical Libraries, Inc., 

Houston, Texas, IMSL Library 2, Reference Manual , pp. FTKALM-1 - 
FTKALM-4, Edition 6, 1977. 



101 



INITIAL DISTRIBUTION LIST 



1. Defense Documentation Center 
Cameron Station 
Alexandria, Virginia 22314 

2. Library, Code 0142 
Naval Postgraduate School 
Monterey, California 93940 

3. Department Chairman, Code 62 
Department of Electrical Engineering 
Naval Postgraduate School 
Monterey, California 93940 

4. Professor H. Titus, Code 62Ts 
Department of Electrical Engineering 
Naval Postgraduate School 
Monterey, California 93940 

5. Robert H. Williams 5 

Naval Weapons Center 

Code 62201 

China Lake, California 93555 



No. Copies 
2 

2 

1 



102 






# 



* 






Thes i s 

w6l8 

c.l 



187022 



Wniiams, 

Two dimensional 
matched filtering for 
attitude measurement 
with applications to 
target tracking. 

l 90 ECeS 27895 



Thesi s 
W6l8 
c.l 



187022 



Willi ams 

Two dimensional 
matched filtering for 
attitude measurement 
with applications to 
target tracking. 



